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Chapter 1 

General Features of CometBoards 




Issues and Strategies in Solving Multidisciplinary 
Optimization Problems 

Surya Patnaik 
Ohio Aerospace Institute 
Brook Park, Ohio 44142 


Abstract 

Optimization research at NASA Glenn Research Center has addressed the design of structures, aircraft 
and airbreathing propulsion engines. The accumulated multidisciplinary design activity is collected 
under a testbed entitled COMETBOARDS. Several issues were encountered during the solution of the 
problems. Four issues and the strategies adapted for their resolution are discussed. This is followed by a 
discussion on analytical methods that is limited to structural design application. 

An optimization process can lead to an inefficient local solution. This deficiency was encountered 
during design of an engine component. The limitation was overcome through an augmentation of 
animation into optimization. Optimum solutions obtained were infeasible for aircraft and airbreathing 
propulsion engine problems. Alleviation of this deficiency required a cascading of multiple algorithms. 
Profile optimization of a beam produced an irregular shape. Engineering intuition restored the regular 
shape for the beam. The solution obtained for a cylindrical shell by a subproblem strategy converged to 
a design that can be difficult to manufacture. Resolution of this issue remains a challenge. The issues 
and resolutions are illustrated through a set of problems: Design of an engine component, Synthesis of a 
subsonic aircraft, Operation optimization of a supersonic engine, Design of a wave-rotor-topping device, 
Profile optimization of a cantilever beam, and Design of a cylindrical shell. This chapter provides a 
cursory account of the issues. Cited references provide detailed discussion on the topics. 

Design of a structure can also be generated by traditional method and the stochastic design concept. 
Merits and limitations of the three methods (traditional method, optimization method and stochastic 
concept) are illustrated. In the traditional method, the constraints are manipulated to obtain the design 
and weight is back calculated. In design optimization, the weight of a structure becomes the merit 
function with constraints imposed on failure modes and an optimization algorithm is used to generate 
the solution. Stochastic design concept accounts for uncertainties in loads, material properties, and other 
parameters and solution is obtained by solving a design optimization problem for a specified reliability. 
Acceptable solutions can be produced by all the three methods. The variation in the weight calculated by 
the methods was found to be modest. Some variation was noticed in designs calculated by the methods. 
The variation may be attributed to structural indeterminacy. It is prudent to develop design by all three 
methods prior to its fabrication. The traditional design method can be improved when the simplified 
sensitivities of the behavior constraint is used. Such sensitivity can reduce design calculations and may 
have a potential to unify the traditional and optimization methods. Weight versus reliability traced out an 
inverted-S-shaped graph. The center of the graph corresponded to mean valued design. A heavy design 
with weight approaching infinity could be produced for a near-zero rate of failure. Weight can be 
reduced to a small value for a most failure-prone design. Probabilistic modeling of load and material 
properties remained a challenge. 
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General Features of CometBoards 

1.0 INTRODUCTION 


Mathematical programming technique of operations research, also referred to as the optimizer, is 
coupled to one or more disciplines to obtain design optimization methodology. It becomes structural 
optimization when the discipline of structures is coupled to an optimizer. Likewise it becomes jet engine 
optimization when the discipline is air-breathing propulsion, and so forth. If the parameters and 
variables of the problem are deterministic in nature then it is referred to as deterministic optimization or 
simply optimization. It becomes stochastic optimization when the parameters become random variables. 
Such variables become distribution functions that are defined by a mean value and a standard deviation. 
Solution to a stochastic design can be calculated for a specified reliability, like one failure in n-number 
of samples. 

A deterministic design- optimization problem can be stated in the following generic from: 

Find n-design variables {X} within upper {X} u and lower{X} L bounds: ({X} 1 " <{X}<{X} U ) 

to optimize an objective function: f(X) 

subjected to m constraints: g y (X)<0 (j = 1, 2, ..., m) 

If there are no constraints then it becomes an unconstrained problem, which has little interest to the 
design-optimization discipline. For stochastic optimization, the deterministic formulation has to be 
modified to account for the uncertainty in load, material properties and sizing variables and the likes. 
There can be variation to statement of the constrained optimization problem given by equation 1 when 
equality constraints are to be added. Design-optimization is discussed first for deterministic case and 
then it is extended into the stochastic domain. 

Optimization research at NASA Glenn Research Center has addressed structural design of 
airbreathing propulsion engines and International Space Station components, aircraft synthesis, as well 
as performance improvement of jet engines. The accumulated multidisciplinary design activity is 
collected under a testbed entitled COMETBOARDS (ref. 1). Several issues were encountered while 
generating solutions to the multidisciplinary problems. This chapter lists four issues and presents the 
strategies that were employed for their resolutions. 

1 . The optimization process produced a local inefficient design for a rear divergent flap of a 
downstream mixing nozzle. An augmentation of animation improved the design. 

2. Solutions obtained for aircraft and engine problems, using single optimization algorithm, 
encountered infeasibility at optimality even though the values of the merit function were in the 
vicinity of the correct solutions. The infeasibility was eliminated through a cascading of multiple 
algorithms. 

3. The shape optimization of a beam produced an irregular shape. The regular shape could be 
restored through engineering intuition. 

4. For a cylindrical shell, the subproblem solution strategy converged to local design that could be 
difficult to manufacture. Resolution of this issue is a challenge. 

This chapter introduces the lessons learned in solving multidisciplinary problems with little emphasis on 
the algorithm or analysis method. The presentation is divided into several topics as listed below. 
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An outline to the optimization testbed COMETBOARDS, 

Stochastic design optimization, 

Description of the illustrative examples, 

Four issues: the local solution and animation, infeasibility and the cascade strategy, 

Irregular design and intuition, substructure solution and manufacturability, multiple design solutions, 
Convergence of optimization algorithms, 

Discussion on analytical design methods, 

Traditional design, Engineering method, Proration technique, 

Discussion on deterministic optimization, 

Stochastic design optimization, 

Comparison of traditional, stochastic and deterministic design methods, 

Conclusions 


2.0 OPTIMIZATION TESTBED COMETBOARDS 

The design optimization testbed COMETBOARDS can evaluate the performance of different 
optimization algorithms and analysis methods while solving a problem. It is a research testbed but not a 
commercial code. The acronym COMETBOARDS stands for “COMparative Evaluation TestBed 
of Optimization and Analysis Routines for the Design of Structures.” The scope of the testbed has 
been expanded to include the design of structures, the synthesis of aircraft, the operation optimization 
of airbreathing propulsion engines as well as stochastic design. COMETBOARDS has three different 
analysis methods and one dozen optimization algorithms. It has a modular organization with a soft 
coupling feature, which allows quick integration of new or user-supplied analyzers and optimizers 
without any change to the source code. The COMETBOARDS code reads information from data fdes; 
formulates design as a sequence of subproblems, and generates the optimum solution. 

COMETBOARDS can be used to solve a large problem, definable through multiple disciplines, each of 
which can be further broken down into subproblems. Alternatively, it can improve an existing system by 
optimizing a small portion of a large problem. Other unique features of COMETBOARDS include 
design variable formulation, constraint formulation, subproblem strategy, global scaling technique, 
analysis approximation through neural network and linear regression method, use of sequential 
or parallel computational platforms, and so forth. The special features and unique strengths of 
COMETBOARDS assist convergence and reduce the amount of CPU time required to solve 
difficult optimization problems of the aerospace industry. COMETBOARDS has been successfully 
used to solve the structural design of the International Space Station components, the design of the 
nozzle components of an airbreathing engine, and airframe and engine synthesis for subsonic and 
supersonic aircraft, mixed flow turbofan engine, wave-rotor-topped engine, and so forth. The 
modular organization of COMETBOARDS is depicted in figure 1 . Brief descriptions of some of its 
modules follow. 
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Scaling and constraint formulation: A multidisciplinary design problem can have a distorted 
design space because its variables and constraints can vary over a wide range. For example, an 
engine thrust design variable, which is measured in kilopounds, is immensely different from its 
bypass ratio, which is a small number. Likewise, the landing velocity of an aircraft measured in 
knots and landing or takeoff field lengths measured in units of thousands of feet differ both in 
magnitude and in units of measure. This module provides a scheme to reduce the distortion by 
scaling the design variables, the objective function, and the constraints such that their relative 
magnitudes during optimization calculations are around unity. The constraints are reformulated to 
alleviate redundancy without affecting the problem definition. Constraint formulation alleviates 
redundancy and reduces their number (ref. 2). The cascade algorithm employs more than one 
optimizer to solve a complex design problem when individual mathematical programming methods 
experience difficulty (ref. 3). 

The module “Analyzers — Structure, Aircraft, Engine — Neural networks and Regression 
approximations” houses three different types of analysis methods. For structural analysis the 
methods available are: COSMIC NASTRAN (ref. 4), MSC/NASTRAN (ref. 5), MHOST (ref. 6), 
Analyze/Danalyze codes (ref. 7), and IFM/ Analyzers (ref. 8). Aircraft analysis can use the FLOPS 
code (ref. 9). The NEPP code (ref. 10) is employed for airbreathing propulsion engine cycle analysis. 
Neural network (ref. 1 1) and regression techniques (ref. 12) can be employed for analysis 
approximation. A problem can utilize any one of three analyzers: (1) an original analyzer, for 
example FLOPS code, or one of the two derived analyzers based on (2) a neural network or (3) a 
regression approximation. 

The module “Engine operations” in figure 1 refers to the performance optimization of airbreathing 
propulsion engines for multiple operation points (ref. 13). “Aircraft synthesis” refers to 
the airframe and engine integration for subsonic and supersonic aircraft (ref. 3). 

The module “Structural designl Subproblem strategy and Parallel computational 
environment” refers to design of structures through regular optimization or a subproblem strategy. 
This strategy is available in sequential and parallel computational environments (ref. 14). “Multiple 
disciplines” refers to the solution of a problem, which is defined through different disciplines. 
COMETBOARDS can accommodate several disciplines each of which can be further divided into 
subproblems. Subproblem strategy is an attempt to alleviate convergence difficulties that can be 
encountered during the solution of a large optimization problem. In this strategy the large problem 
is replaced by a sequence of overlapping modest subproblems. The solution to the large problem is 
obtained by repeating the solution to the set of subproblems until convergence is achieved. 
Substructure strategy, a key module of COMETBOARDS, is further illustrated through the 
design of a cargo-bay support of the International Space Station. The support structure is fabricated 
out of four plates, a cluster of plates referred to as a box, and five beams, as shown in figure 2. For 
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the purpose of design, the support system was divided into four active and one passive substructures. 
The first substructure was the closed box (FGHKIJ) consisting of five plates. Its finite element 
model was obtained by discretizing it into 72 (QUAD4 ) shell elements. The second substructure 
was a trapezoidal plate (FHEC), and its finite element model had 37 shell elements. The third and 
fourth substructures were triangular plates (GHE and GHD) with 12 finite elements each. The fifth 
substructure contained the five connecting beams BD, BG, AD, AG, and AF. The beam designs 
were not changed. These passive variables did not participate in the design calculations but were 
retained during reanalysis. Two independent finite analysis codes LE_HOST and MSC/NASTRAN 
were used to verify the finite element model of the support structure. The finite element model with 
a total of 133 QUAD4 shell and 20 beam elements was considered adequate for design optimization 
because both analyzers produced an acceptable level of accuracy for stress, displacement, and 
frequency. The substructures were clustered next to obtain a set of subproblems. A subproblem contains 
a substructure along with some or all of its neighboring substructures. The four substructures were 
clustered to obtain the following four subproblems: 

Subproblem 1: substructures 3 and 4 Subproblem 3: substructures 1 and 2 
Subproblem 2: substructures 1 and 4 Subproblem 4: substructures 2 and 3 


Notice the coupling between subproblems: substructure 1 is common to subproblems 2 
and 3. Likewise, substructure 2 is common to subproblems 3 and 4, and so forth. Adequate coupling 
between substructures is required for convergence. Inappropriate coupling can increase the amount 
of computation and/or encounter convergence difficulties. At this time substructure coupling is 
decided intuitively. However, it may be possible to automate substructure coupling through the 
gradient scheme developed in reference 15. The COMETBOARDS testbed can optimize a system 
that can be defined in terms of 100 optimization subproblems (ref. 14). In the module “Problem 
formulation and solution,” information is read from data files, the design is cast as a sequence of 
optimization subproblems, and the solution is obtained. The COMETBOARDS testbed is written in the 
Fortran 77 language except for the neural network algorithm, which is written in the C++ language. The 
testbed is available in the Unix operating system in workstations and Cray computers. 


3.0 STOCHASTIC DESIGN OPTIMIZATION 

CometBoards has been extended into the stochastic domain. The stochastic design optimization 
methodology (SDO) in the code has been developed to design components of an airframe structure that 
can be made of metallic and composite materials. The design is obtained as a function of the risk level, 
or reliability, p. The design method treats uncertainties in load, strength, and material properties as 
distribution functions, which are defined with mean values and standard deviations. A design constraint 
or a failure mode is specified as a function of reliability p. Solution to stochastic optimization yields the 
weight of a structure as a function of reliability p. Optimum weight versus reliability p traced out an 
inverted-S-shaped graph. The center of the inverted-S graph corresponded to 50 percent (p = 0.5) 
probability of success. A heavy design with weight approaching infinity could be produced for a near- 
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zero rate of failure that corresponds to unity for reliability p(orp= 1). Weight can be reduced to a small 
value for the most failure-prone design with a reliability that approaches zero {p = 0). Reliability can be 
changed for different components of an airframe structure. For example, the landing gear can be 
designed for a very high reliability, whereas it can be reduced to a small extent for a raked wingtip. The 
SDO capability is obtained by combining three codes: (1) The MSC/Nastran code was the deterministic 
analysis tool, (2) The fast probabilistic integrator, or the FPI module of the NESSUS software (ref. 16), 
was the probabilistic calculator, and (3) NASA Glenn Research Center’s optimization testbed 
CometBoards became the optimizer. The SDO capability requires a finite element structural model, a 
material model, a load model, and a design model. The stochastic optimization concept is illustrated 
considering an academic example and a real-life raked wingtip structure of the Boeing 767-400 
extended range airliner made of metallic and composite materials. 


Optimization Algorithms in COMETBOARDS: One dozen mathematical programming algorithms and a 
cascade strategy are available in CometBoards. A list of the algorithms in no particular order follows. 

(1) Method of feasible direction (FD) (ref. 17) 

(2) Modified feasible direction method (mFD) (ref. 1 8) 

(3) Reduced gradient method (RG) (ref. 19) 

(4) Generalized reduced gradient method (GRG) (ref. 20) 

(5) Sequence of unconstrained minimization technique (SUMT) (ref. 21) 

(6) Sequential linear programming (SLP) (ref. 17) 

(7) Sequential quadratic programming method (SQP) (ref. 22) 

(8) IMSL optimization routine (IMSL) (ref. 23) 

(9) NPSOL package of NAG library (NLPQ) (ref. 24) 

(10) Sixteen different versions of optimality criteria methods (OC) (ref. 25) 

(1 1) A genetic algorithm (GENMO) (ref. 26) 

(12) Fully utilized design algorithm (FUD) (ref. 15) 

Different algorithms employ different strategies to calculate the search direction and the step 
length. The comparative performance of the algorithms in solving a set of 45 problems is reported 
in reference 27. The performance of six algorithms in solving a set of medium and large structural 
problems is depicted in a bar chart in figures 3(a) and 3(b), respectively. Examples with design 
variables in the range of 20 to 39 with about 200 behavior constraints are referred to as medium 
problems, see figure 3(a). Examples with more than 40 independent design variables and several 
hundred constraints are referred to as large problems, see figure 3(b). The success rate of different 
optimization algorithms for 10 large structural problems is also depicted in a Venn diagram in 
figure 3(c). Cascade solutions for the problem are given in a table (see insert 3(e)). The success 
of an algorithm is represented by unity, which is the normalized value of the merit function, see 
figures 3(a) and 3(b). A normalized value of less than unity (an infeasible solution) or greater than 
unity (an overdesign condition) represents underperformance. Most optimizers available in the testbed 
solved at least one -third of the examples, but none of the optimizers could successfully solve all 
the problems. Every structural problem could be solved by at least one of the six different optimization 
algorithms. However, even the most robust optimizer encountered difficulty in generating optimum 
solutions for aircraft and engine problems. 
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4.0 ILLUSTRATIVE EXAMPLES 


For the purpose of illustration, we have selected six problems: 

Problem 1 — Structural design of an engine component. 

Problem 2 — Synthesis of a subsonic aircraft. 

Problem 3 — Operation optimization of a supersonic engine. 

Problem 4 — Design of a wave-rotor-topping device. 

Problem 5 — Profile optimization of a cantilever beam. 

Problem 6 — Design of a cylindrical shell. 

Augmentation of animation into optimization is illustrated through Problem 1. Problems 2, 3, and 4 
illustrate the cascade strategy. Problem 5 addresses shape optimization. The subproblem solution 
strategy is illustrated through problem 6. Brief definitions of the problems follow. 

4.1 Problem 1: Structural Design of an Engine Component 

A mixed-flow turbofan engine exhaust nozzle referred to as a “downstream mixing nozzle” 
for a High Speed Civil Transport aircraft to operate at a cruise speed of Mach 2.4 and in a range of 
5000 nautical miles is shown in figure 4(a). It is fabricated out of rear and forward divergent flaps, 
rear and forward sidewalls, bulkheads, duct extensions, six disk supports, and other components. 

The design complexity of the nozzle is increased with flight Mach number, pressure ratio, temperature 
gradient, dynamic response, and degradation of material properties at elevated temperature. 

The flap is made of Rene 125 material with a Young’s modulus of 30.4 million psi, a Poisson 
ratio of 0.3, a density of 0.308 lb/in.3 and an allowable strength of 1 17 ksi. The flap shown in 
Figure 4(b) has a length of 96 in. and a width of 72 in. It is supported by two variable-depth edge 
beams with a maximum depth of 14 in. A grid with a spacing of 12 in. by 12 in. and a depth of 2 in. 
stiffens the flap. MHOST and MSC/NASTRAN analyzers were used for the static as well as the 
dynamic analyses of the flap. A finite element model with 5593 degrees of freedom and 946 
QUAD4 elements was used for analysis because both methods produced acceptable values for 
stress, deformation, and frequency. For this problem, the thickness of the edge beams, stiffeners, 
and skin panel were considered as the design variables. The flap was designed for minimum weight 
for von Mises stress, local buckling, displacement, and frequency constraints. The problem had 946 
stress and 946 instability constraints, as well as four displacements and one-frequency constraints. 

4.2 Problem 2: Synthesis of a Subsonic Aircraft 

The system synthesis capability of COMETBOARDS is illustrated through a subsonic 
aircraft to operate at Mach 0.85 cruise speed. Solution of this problem required a soft coupling of 
COMETBOARDS to NASA Langley’s aircraft analysis code FLOPS. The FLOPS analyzer 
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encompasses several disciplines: weight estimation, aerodynamic analysis, engine cycle analysis, 
propulsion data interpolation, mission performance, airfield length requirement, noise footprint 
calculation, and cost estimation (refs. 28, 29, and 30). The objective of the synthesis problem was 
to determine an optimum airframe-engine design combination for a set of active design variables 
under specified engine and aircraft behavior parameters to minimize a composite merit function that 
could be generated as a linear combination of weight, mission fuel, lift-to-drag ratio, and NOx 
emission. The design variables considered were engine thrust, wing area, engine design point 
turbine entry temperature, overall pressure ratio, bypass ratio, and fan pressure ratio. Constraints 
were imposed on mixed approach climb thrust, second segment climb thrust, landing approach 
velocity, takeoff field length, jet velocity, and compressor discharge temperature. For the subsonic 
aircraft model, the gross takeoff weight was considered as the merit function. 

4.3 Problem 3: Operation Optimization of a Supersonic Engine 

The operation optimization of airbreathing propulsion supersonic engine required a soft 
coupling of the analyzer NEPP, the NASA Engine Performance Program, to the optimization 
testbed COMETBOARDS. The engine-cycle NEPP code can configure and analyze almost any 
type of gas turbine engine that can be generated through the interconnection of a set of standard 
physical components such as propellers, inlets, combustors, compressors, turbines, heat exchangers, 
flow splitters, nozzles, and others. An engine can be designed for different types of hydrocarbon jet 
fuels, cryogenic fuel, and slurries. For their thermodynamic analysis, built-in curve fits generated 
from empirical results have been incorporated into the NEPP code. A description of the NEPP code 
with typical input files for a set of engine configurations can be found in references 3 1 and 32. 

The engine operation optimization problem, with its associated design variables (such as the 
engine shaft speed, the wave-rotor speed, the flow area, the geometrical parameters of the ducts, 
etc.) and constraints (imposed on pressure ratios, surge margins, temperature limits, and entrance 
Mach number, etc.) was cast as a sequence of nonlinear optimization subproblems with thrust as the 
merit function. Engine operation design became a sequence of interdependent problems, or one 
optimization subproblem for each operating point. The optimization process typically adjusted a 
few engine parameters. The difficulty in the engine problem did not lie with the number of active 
design variables, but it was associated with its multiple operating-point character, constraint validity 
ranges, and the iterative nature of engine cycle analysis. The most reliable individual optimization 
algorithm available in COMETBOARDS could not produce a satisfactory feasible optimum solution 
for this engine problem because of the large number of operation points in the flight envelope, 
the diversity of the constraint types, and the overall distortion of the design space. However, 
COMETBOARDS’ unique featuresl which included a cascade strategy, variable and constraint 
formulations, and scaling devised especially for difficult multidisciplinary applications 
successfully optimized the performance of subsonic and supersonic engine concepts. Even when 
started from different design points, the combined COMETBOARDS and NEPP codes converged to 
about the same global optimum solution. This reliable and robust design tool eliminated manual 
intervention in the design of the airbreathing propulsion engines and eased the cycle analysis 
procedures. Engine design is illustrated through a supersonic Mixed-Flow Turbofan Engine (MFTF) 
problem. It was configured with 15 components and was designed for a flight envelope with 122 
operating points. The design required the solution of a sequence of 122 optimization subproblems. 
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The objective was to maximize the net thrust of the supersonic engine at each operating point, 
accounting for an installation drag. Each subproblem had 3 independent design variables and 22 
behavior constraints. 

4.4 Problem 4. Design of a Wave-Rotor-Topping Device 

A high bypass-ratio subsonic wave-rotor topped turbofan engine is made of 16 components 
that are mounted on two shafts with 21 flow stations. It was modeled with standard components 
that include an inlet and a splitter, then branching off to a compressor, a duct, and a nozzle. The 
main flow proceeded through a fan, a duct, a high-pressure compressor, a duct, a high-pressure 
turbine, a low-pressure turbine, a duct, and a nozzle. The components mounted on the first shaft 
included the fan, the compressor along the secondary flow branch, the low-pressure turbine, and a 
load. The second shaft carried the high-pressure compressor along the main flow, the high-pressure 
turbine and a load. The four-port wave-rotor (with burner inlet and exhaust, compressor inlet, and 
turbine exhaust), was located between the high-pressure compressor and the high-pressure turbine. 

The engine operating condition was specified by a 47-mission flight envelope, with altitude in the 
range (sea level to 40 000 ft) and speed between (0.0 to 0.85) Mach. To examine the benefits that 
accrue from the wave-rotor device, the engine was optimized considering several of its baseline 
variables and constraints, not explicitly associated with the wave-rotor, as passive. The design 
objective was to maximize the net engine thrust at each of the 47 operating points. It had two 
design variables: heat added to the wave-rotor and the wave-rotor speed. Its 16 behavior constraints 
included the corrected speed ratio for the compressor and the fan, the unmixed wave-rotor temperature, 
the surge margin on the compressor, and the fan pressure ratio for the turbine. 

4.5 Problem 5: Cantilever Beam 

Calculation of an optimum profile or depth along the length is illustrated considering the 
cantilever beam shown in figure 5. The beam is made of aluminum with a Young’s modulus 

E = 10xl0 6 ps/), a Poisson’s ratio (t> = 0.3) , and a weight density (p = 0.1/6 / w 3 ). It is 240 in. long 

and 6 in. deep. The beam weight was the merit function and stress and displacement were its behavior 
constraints. For the purpose of analysis, the beam was modeled with six 20-node-hexahedral 
elements of the IFM/ Analyzers code (ref. 33). The finite element model with 160 displacement 
degrees-of-ffeedom was used for the analysis model because it adequately predicted the stress and 
displacement responses of the beam. The profile optimization problem had nine depth design 
variables, and its nine constraints included eight stress and one displacement limitations. 


NASA/CR— 201 3-2 1 7748 


11 



4.6 Problem 6: Cylindrical Shell 


A cylindrical shell with rigid diagrams under two line loads is shown in figure 6. It is made of 
steel with a Young’s modulus {E = 30x1 0 6 psi } , a Poisson’s ratio (u = 0.3) , and a weight density 

[p = 0.289 lb I m 3 )lb/in3. It has a radius of 100 in. and length of 200 in. Because of symmetry, only one- 

eighth of the structure was considered for design. A finite element model with 100 QUAD4 elements of 
the MHOST analyzer was considered adequate to predict its response. Its design was cast as an 
optimization problem with weight as the objective function and constraints were imposed on stress and 
displacement. The thickness of the cylinder along its length and circumference were considered as 
design variables through a profiled depth link factor to provide a heavier design under the load. The one 
eighth shell model was divided into four substructures along its length. The substructures were clustered 
to obtain three and four subproblems for sequential and parallel computations, respectively. An 
additional fourth subproblem was required to avoid convergence difficulty in parallel calculations. 


5.0 Issues in Multidisciplinary Design Optimization 

The four issues: (a) local solution and animation, (b) infeasibility and cascade strategy, (c) irregular 
design and intuition, and (d) substructure solution and manufacturability are addressed in this section. 

5.1 Local Solution and Animation 

The design of the flap or Problem 1, which is depicted in figure 4(b), was obtained using 

three optimization algorithms. All three methods produced the same optimum weight of 1448.2 lb. 

At the optimum, the frequency at 40 Hz and displacement at 1 in. were the active constraints. Stress 
and local instability were passive constraints. Its animation at the 40 Hz frequency was examined, 
and one frame is depicted in figure 4(c). The animation exhibited a local frequency condition. The 
edge beams vibrated with significant amplitude, while the response of the rest of the structure was 
insignificant. In other words, only a small part of the structure carried a major portion of the load. 

In this design, the edge-beams became more failure prone than other parts of the structure. 

The configuration of the flap was modified through an examination of the animation of a 
series of designs. The final modified configuration was obtained by increasing the depth of a single 
edge stiffener to 6 in. from its original 3 -in. depth. This configuration was optimized. The increase 
in the material for this stiffener was compensated for by a reduction in the thickness of the edge 
beam by 58 percent between the two configurations. The dynamic animation of the flap and the von 
Mises stress distribution are shown in Figures 4(d) and 4(e), respectively. The structure vibrated in 
a breathing type of mode, or the entire flap responded as a single unit. The optimum design also 
displayed a full stress condition for a major portion of the flap as shown in figure 4(e). The optimum 
weight at 1204.7 lb. was 20 percent lighter than the original configuration. Animation assisted 
optimization reduced the weight and improved the overall design of the flap. 
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5.2 Infeasibility and Cascade Strategy 


Individual optimization algorithms encountered difficulty in generating solutions to aircraft 
and engine problems. A useful strategy that combined the strength of more than one optimizer 
was conceived. This cascade strategy, using a number of optimization algorithms, one followed 
by another in a specified sequence, was found to be superior to the component optimizers, see 
figure 3(d). The cascade algorithm was employed to solve the 10 large structural problems referenced 
in Section II. 1 . The success rate of the cascade strategy and the individual algorithms is 
shown in Figures 3(c) and 3(e). The cascade strategy solved all 10 problems. 

Cascade solutions to the subsonic aircraft, the supersonic mixed-flow turbofan engine, and 
the subsonic wave-rotor- topped engine are shown in figure 7. The subsonic aircraft system design 
problem was solved using a three-optimizer cascade algorithm: optimizer 1, followed by optimizer 2 
and optimizer 1 again. The cascade solution, along with solutions obtained from individual algorithms, 
is depicted in figure 7. Optimizer 1 , when used alone, converged to a heavy infeasible solution 
at 202 005 lb for the aircraft weight see figure 7(a). Likewise, optimizer 2 alone also produced a 
heavy design at 202 854 lb, see Figure 7(b). Optimizer 1 required a takeoff field length of 6282 ft 
against an available 6000-ft runway. Optimizer 2 overestimated the excess fuel requirement at 8150 gal 
against a desired amount of 5000 gallons. The convergence rate of the two algorithms differed 
producing solutions that were 1 percent and 2 percent overweight, respectively, for the two 
algorithms, see figures 7(a) and 7(b). The two solutions, although close to the optimum, were not 
attractive to industry because they violated the safety margins. A cascade algorithm created from 
the same two optimizers (optimizer 1 - optimizer 2 - optimizer 1) successfully solved the problem 
with a feasible optimum solution at 199 818 lb for the aircraft weight, see figure 7(c). 

Solution of the mixed- flow supersonic and wave-rotor-topped subsonic engines also 
required the cascade strategy. A two-optimizer cascade, (optimizer 3 followed by optimizer 2) 
successfully solved the supersonic engine problem, see figure 7(d). The solution to the wave- 
rotorenhanced subsonic engine, see figure 7(e), required a three-optimizer cascade algorithm: optimizer 
1 followed by optimizer 2 and optimizer 1 again. Aircraft (FLOPS) and engine (NEPP) analyzers are 
nonlinear codes that incorporate multiple disciplines. The analysis assumptions, dependent on altitude, 
Mach number, and engine power setting, can be challenged. The Newton-Raphson iterations during 
“engine-balancing” may not converge, leading to an engine “stall condition” with zero thrust, or to the 
NEPP code producing a NaN (not a number) for some constraints. A single optimization algorithm can 
terminate abruptly, or hit a constraint boundary without any improvement to the merit function. To 
overcome the difficulties in aircraft and engine analysis codes we employed two competing 
approximators: linear regression and neural network methods. Even then, the cascade algorithm was 
required because individual algorithms encountered difficulty (ref. 13). The nature of engine and aircraft 
problems required the strength of multiple algorithms or the cascade strategy even with the analysis 
approximators, which however, was not expected. 


5.3 Irregular Design and Intuition 

The profile optimization of the cantilever beam, or Problem 5, was attempted with uniform 
upper and lower bounds for the design variables at 0.5 in. and 15 in., respectively. The 
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optimization problem was solved using optimizer 4. The solution obtained is shown in figure 5. 

The optimum weight was 1697.5 lb, and the root stress was the only active constraint. An 
odd-shaped profile shown in figure 5 was obtained. The profile was peculiar because the free end, 
corresponding to a zero stress condition, had a depth of 5.14 in. instead of an anticipated lower 
bound depth of 0.5 in. The optimum solution most likely challenged linear structural analysis 
assumptions, and the optimization iterations encountered difficulties. The situation did not improve 
when a different optimization algorithm or a cascade strategy was employed. The problem was 
solved successfully with manual intervention. Instead of uniform upper or lower bounds and a 
single optimization algorithm, the design was cast as a sequence of subproblems. Upper and lower 
bounds and initial designl through engineering intuitionl were progressively changed for the 
subproblems. This procedure produced a converged solution that is shown in figure 5. It had the same 
minimum weight of 1697.5 lb, which was identical to the odd-shaped design. Its root stress and tip 
displacement were the active constraints. For this problem, optimization algorithms converged to 
two distinct local solutions with equal minimum weight. Industry will be more inclined to accept 
the monotonically profiled beam than the odd-shaped design. 

The flap design and beam profile calculations represent typical problems of the aerospace 
industry. Neither problem could be solved satisfactorily without manual intervention. The difficulties 
encountered to some extent justify the reluctance of the aerospace industry to accept 
advanced optimization methods, abandoning their time-tested design rules. Neither mathematical 
programming algorithms nor the traditional design rules could produce optimum hardware 
configurations that could be manufactured. Their combination, however, was a winner. 

5.4 Subproblem Solution and Manufacturability 

A no-manufacturable solution obtained in the subproblem strategy is illustrated considering 

the design of a cylindrical shell problem, or Problem 6. Solution to the problem was obtained using: 

(1) Regular optimization, where the entire structure was considered as a single problem. 

(2) Subproblem solution, wherein four overlapping substructures were used. 

The optimum profiles for the cylindrical shell obtained using regular and subproblem 
strategies are depicted in figure 6. The two optimum weights obtained were 1 161.95 lb and 1 154.1 
lb for regular and subproblem strategies, respectively. The difference of 0.676 percent in the 
solutions could be considered negligible especially for the problem complexity. The depth differed 
substantially between the two solutions. At the crown, the optimum depths of 1.322 in. and 2.471 
in. obtained by the two methods varied by 53 percent. At the optimum, the regular optimization and 
the subproblem strategy produced a different number of active constraints. The subproblem strategy 
produced only active stress constraints, whereas both stress and displacement constraints were 
active for the regular optimization. The profde obtained using regular optimization was more uniform 
than that generated through subproblem optimization. To verify the existance of two different 
designs (one obtained from the subproblem strategy and the other from the regular optimization), we 
resolved the regular optimization case by setting the initial design equal to the optimum solution 
that was generated from the subproblem strategy. The solution converged to the initial design, 
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confirming the existence of the two local solutions with about the same minimum weight. For this 
problem the attractiveness of subproblem strategy has been challenged because the design thus 
obtained is more difficult to manufacture compared to the regular solution. The authors are not 
aware of a scheme to alleviate this limitation of the subproblem solution strategy. 

5.5 Multiple Design Solutions 

More than one optimum solution can be obtained for engineering design problems. For 
example, different solutions were obtained for the nozzle flap, beam profile, cylindrical shell, subsonic 
aircraft, and engine problems. The values of merit function changed little between the different 
solutions. The intermediate designs, however, could become nonmanufacturable or their safety could be 
questioned. These limitations were alleviated through animation and engineering intuition. The 
utilization of animation improved the design of the nozzle flap, and manipulation of the design bounds 
promoted the manufacturability of the beam. Philosophically, we can arrange design methods into a 
spectrum. Experience occupies the bottom strata. Optimization methods belong to the upper spectrum. 
Popular design rules occupy the central spectrum. Design improved as we moved from the lower to 
upper spectrum methods but complexity increased. The design optimization method is in existence for 
about half a century, yet its full potential has yet to be exploited by industry. The situation can be 
improved by merging design optimization with the engineering knowledge contained in the broad 
spectrum methods. This goal can be achieved when proponents of the optimization method and 
industrial designers work in tandem. 

5.6 Convergence of Optimization Algorithms 

Consider the convergence of the subsonic aircraft weight by two different algorithms, shown in 
Figures 7(a) and 7(b). Convergence was monotonic by both methods even though the rate differed, as 
expected. A similar trend was observed for the engine problems. Convergence at times oscillated 
about a mean solution until the maximum iteration limit was reached. Redundancy of the active 
constraints, their continuity and convexity may have created this situation. Engineering design problems 
may not satisfy some of the basic requirements that form the foundation of mathematical programming 
methods. For example, the stress and displacement constraints of structural problems are by nature 
redundant (ref. 2). The specified range of the constraints of engine and aircraft problems are 
susceptible to violation during optimization calculations. Optimum solution to multidisciplinary 
problems of the industry have to be obtained utilizing available analysis and optimization capabilities. 
The cascade strategy was found successful in generating reliable solutions for structures, aircraft, and 
engine problems. The cascade strategy was required even when neural network and regression 
approximators were used to approximate constraints and merit functions of the engineering design 
problems this however was not expected. 

An animation-assisted optimization successfully solved the flap design problem. Cascading of multiple 
algorithms solved aircraft and engine problems. The beam profile problem was solved through an 
incorporation of engineering intuition. Generating a regular manufacturable design using a subproblem 
optimization scheme still remains a challenge. Bringing optimization methods to their rightful 
industrial environment from the academic plane requires the combined effort of designers and 
optimization researchers. 
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6.0 ANALYTICAL DESIGN METHODS 


Analytical methods to design a structure can be grouped under three categories: 

1. Traditional design 

2. Optimization method and 

3. The stochastic design concept 

The traditional method practiced in the industry is based on fully stressed and fully utilized concepts. 
The stress and stiffness constraints are manipulated until a satisfactory design is obtained. The weight of 
the structure is not used explicitly in traditional design calculations but it is back-calculated. A design is 
generated in few iterations and gradient or sensitivity of weight or constraint is not used. 

In design optimization method the weight of a structure becomes the merit function and constraints are 
imposed on failure modes, etc. An optimization algorithm of mathematical programming method of 
operations research is used to generate the solution to the design problem. The optimization method is 
based on advanced calculus and it is highly mathematical in nature. 

Stochastic design concept accounts for uncertainties in loads, material properties, design variables and 
failure modes. Probabilistic distribution functions are used to accommodate the uncertainties and 
response is estimated for a specified reliability. Design is generated for the specified reliability, like one 
failure in n-number of samples and an optimization algorithm is used. Modeling of parameters of 
structural design as distribution functions can be a challenge. This method is computation intensive. 

A preliminary comparison of the three methods is attempted in this section. Similarities and differences 
are examined. Numerical solutions included an industrial strength design problem in addition to 
academic examples. Solution to numerical examples utilized three codes. MSC/Nastran code was the 
deterministic analyzer. Fast probability integrator or the FPI module of the NESSUS software was the 
probabilistic response calculator. NASA Glenn Research Center’s testbed CometBoards became the 
optimization tool. 


6.1 Formulation of Structural Design Methods 

The three analytical design methods namely, the traditional design, optimization method and the 
stochastic design concept were formulated considering the example of a tapered beam as an illustration. 
The beam shown in figure 8 was made of aluminum material. The determinate structure was subjected 
to a tip load of 20 kip. The beam depths (di and di) were the two design variables while the passive 
thickness variable (f = 4.0 in.) was uniform across its length. Other parameters were shown in the figure. 
The design objective was to minimize weight under strength and stiffness constraints. For stochastic 
optimization, normal distribution was considered for the primitive random variable with a mean value 
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(//) and a standard deviation (cr) . It was assumed that mean value (//) was close to the nominal value 
that was used in deterministic design calculation. Standard deviation O) was considered a percent of 
mean value. 


6.1.1 Traditional Design Method 

The traditional design method, practiced in the industry, is based on a fully utilized concept (ref. 34). 
Design is obtained for strength limitation and it is then adjusted for violated stiffness constraints, if any. 
For the determinate beam in figure 8 the two depths were calculated from the induced stresses at 
locations 1 and 2 as follows: 


'mA 

W Ji 


and 



(i) 


where, a 0 was the stress allowable, induced stress was a, the bending moment was M, moment of inertia 
was 7, and y represented the half-depth of the beam. Location 1 is at the clamped boundary, while 
location 2 is at the transition, see figure 8. 


Solution to equation 1 yielded the traditional design for the determinate beam. Its weight was back 
calculated and it could be close to the minimum weight even though it was not explicitly used in the 
design calculation. For an indeterminate structure, the equation set 1 have to be solved iteratively. Even 
for a determinate structure, the design has to be updated when a displacement limitation is violated. It is 
not straightforward to redesign for violated displacement constraint (ref. 15). Two approaches are 
discussed. These are the engineering method and the proration technique. 


6.1.2 Engineering Method 

The experience based ‘Engineering Method’ attempts to satisfy the violated displacement constraint 
intuitively. For the tapered beam, the displacement (8) at the free end can be expressed in closed form 
as: 
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In the engineering method, all depths are changed uniformly by updating the depths by a factor (k) as: 


d" ew = k d* 
d n 2 ew = k d s 2 


( 3 ) 
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The design for strength constraints only is designated with the superscript ‘s’ in equation 3. The design 
with superscript ‘new’ satisfies both strength and stiffness constraints, assuming k to be greater than 
unity (k > 1.0) . Solution of equations (2 and 3) for the required displacement limit (£ max ) yields the 


factor (k) as: 
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Equation 4 can be modified and solved iteratively for a situation with multiple violated displacement 
constraints. Pick the most violated displacement limitation and make it feasible using the procedure 
discussed. Reevaluate the design to determine its feasibility. Satisfy a next violated displacement 
limitation if any. A feasible design is generated in a few iterations. For a general type of structure, a 
solution to equation 4 or its variant is quite simple in a computer. Engineering method may exhibit a 
tendency to become heavy because the design update is uniform across all members. A satisfactory 
design for the problem may be obtained by updating the root depth (c/J without any change to (c/ 2 ) . 

This concept of selective updating of the design is discussed next. 


6.1.3 Proration Technique 


The proration technique is an attempt to alleviate the limitation of the engineering method, which may 
have a tendency to produce an over-design condition. Instead of uniform proration, a weighted 
proration can be used. The weight factors are calculated from the sensitivity of the displacement 
constraint. For the tapered beam the sensitivity can be calculated in closed form as: 
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{0.36, normalized to pr : = 1 .75} 
{0.19, normalized to pr 2 = 1 .00} 


( 5 ) 


Then equation 3 can be rewritten as: 


d \ new = (pp = 1 .75 )k cf 
d n 2 ew = (pr 2 =1.00 )kd s 2 


The weight factors bias the design towards (d , ) which is the depth of the cantilever near the support and 

it is intuitively correct. The proration method is discussed in reference 2 for indeterminate trusses. The 
calculation of sensitivity has been simplified and it does not require intensive computation (ref. 35). 
The engineering method and the proration technique were developed considering the example of a 
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determinate structure. The formulations, however, can be used for indeterminate flexural structures. 
Design iterations are required for strength and stiffness limited design for an indeterminate structure. 
The infeasibility of a fully stressed design (ref. 2) has to be addressed when the formulation is applied to 
indeterminate trusses. This design deficiency may not be pronounced for indeterminate flexural 
structure. 

6.2 Optimization Method 

In this popular method structural design is cast as a mathematical programming problem with weight as 
the merit function and failure modes (here strength and stiffness limitations) as the constraints. Earlier, 
design was cast as an optimization problem. The design optimization problem is solved using 
mathematical programming techniques of operations research. We solve the problem utilizing the 
optimization test-bed developed at NASA Glenn Research Center. The code, referred to as 
CometBoards formulates structural design and solves it as a nonlinear mathematical programming 
problem. Problem solution can use any one of a dozen optimization algorithms or the cascade strategy. 
CometBoards has been successfully applied to design structures, of air-breathing propulsion engines, 
airliner synthesis and other problems. 

6.3 Stochastic Design Optimization 

In stochastic design optimization the variables, merit function and the constraints are treated as 
random parameters. Such a parameter is defined through a distribution function, (like normal, Weibull, 
lognormal, and Gumbel etc.), with a mean value and a standard deviation. The values of the random 
parameters can be estimated for a specified reliability (p) like one failure in n number of samples. The 
design variables, merit function and constraints become function of the reliability p. In stochastic 
optimization, a design became a function of risk, or reliability, p. The design method accommodates 
uncertainties in load, material properties, and failure theory and design variables. Solution of stochastic 
optimization yields the design and weight of a structure as a function of reliability p. When optimum 
weight versus reliability p is plotted, an inverted-S-shaped graph is obtained. The center of the inverted- 
S graph corresponds to a 50 percent probability of success. A heavy design with weight approaching 
infinity could be produced for a near-zero rate of failure that corresponds to unity for reliability {p = 1). 
Weight can be reduced to a small value for the most failure-prone design with a reliability that 
approaches zero {p = 0). 

In stochastic optimization (ref. 36), a constraint gj for a specified percent probability (p between 
0 and 1), is generated from random response variables within prescribed random upper and lower 

bounds g U j and gj, respectively, as: (gj < gj < g^ )> p. A typical stochastic constraint (for stress 

limitation r/ < r/o) has the following form: 
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where, p is reliability, p n is the mean value of stress, p is the allowable stress, 0 * is the critical phi 

parameter, cr is the standard deviation of stress and a TO is the standard deviation of stress allowable. A 

stochastic optimization capability was developed by combining three codes: The MSC/Nastran code was 
the deterministic analysis tool, the fast probability integration code, or the FPI module of the NESSUS 
software was the probabilistic calculator and NASA Glenn Research Center’s optimization testbed 
CometBoards became the optimization tool. 


6.4 Comparisons of the Design Methods 

Several examples were solved by traditional, optimization and stochastic design methods. The design 
solutions were compared for weight, design variables and constraints. Three examples discussed are: a 
determinate tapered beam, an indeterminate beam and the raked wingtip of the Boeing 767-400 
extended range airliner made of metallic and composite materials. Solution used the MSC/Nastran code 
as the deterministic analyzer. Fast probability integration or the FPI module of the NESSUS software 
was the probabilistic response calculator. NASA Glenn Research Center’s testbed CometBoards was 
the optimizer. 

6.4.1 Tapered Beam 

The tapered beam shown in figure 8 was made of aluminum with a nominal value for Young’s modulus 

3 

of(E = 10 ksi ) , a weight density (p = 0.1 lb I in ) and a length of 2 a = 144/n. The uniform thickness of 
the beam was set to (f = 4 in). The beam was subjected to a transverse load ( P = 20 kip) at the tip. The 
allowable stress limit was (<r 0 = 20 ksi) and the tip displacement was not to exceed one inch (s < 1) . The 
depths (d v d 2 ) were the two design variables. Minimum weight design was the objective. The random 

variables for stochastic calculations were listed in table 1 . Normal distribution was assumed for all 
random variables. Their mean values were set equal to or close to the nominal values with some 
variations. For example, for the nominal strength at 20 ksi, a mean value was set at 25 ksi with a 
variation of 2.5 ksi. Likewise, other random parameters were defined in table 1. Specification of realistic 
mean values and deviations was a challenge. Those in table 1 were assumed for subsequent calculations. 
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Table 1. Mean value and standard deviation for random variables. 


Primitive variable 

Mean value 

Standard deviation 

Stress limit: <x n 

u = 25 ksi 

'< T n 

<7 = 2.5 ksi 

a o 

Disp. Limit: S 

p s = 1 .25 in 

u s =0.125 in. 

Modulus: E 

p E = 10 ksi 

o E =\ ksi 

Density: p 

A, =0.1 

o D = 0.005 in. 

Load: P 

p P =15 kip 

<j p = 3 kip. 

Depth: d, and d. 

/u d = 25 in and ju d2 =10 in 

<r di = 3 in and o d2 = 2 in 

Thickness t 

''si- 

ll 

cr t = 0.25 in.. 


Both deterministic and stochastic designs were generated following the methodology discussed 
earlier. For the strength-limited design, the traditional stress ratio method yielded the two depths, that 
were calculated from moments at the base and at the transition point for the nominal strength of 
(<r 0 = 20 ksi). The optimization methods also converge to the same design point. The design 

(d, = 14.7/n., d 2 = 1 0.4 in. and weight = 722.9 Ibf ) with a tip displacement of (d max = 2.308 in.) violated the 
stiffness limitation of one inch(<5 < 1) . 


The design had to be adjusted for the feasibility of the displacement limitation. Generation of a 
design for both stress and stiffness constraints was not straight forward even for the simple determinate 
structure. The maximum displacement (<5 max ) , which occurs at the tip of the cantilever, can be expressed 

in terms of the two depths as: 


: 746.5 
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The displacement function was graphed in figure 9. The x-axis represented a pair of (d, and d 2 ) or the 
design for which (d max = 1 in) . It was observed that there were an infinite number of designs that satisfied 
the displacement limitation. Designs generated by different methods are marked in figure 9. 

For the problem, the optimality criteria method produced a design with the least weight of 946.9 lbf. 
The weight calculated by the prorated design method at 955.4 lbf was about one percent heavier than the 
solution by the optimality criteria approach. The design was heavier by 7.5 percent by the engineering 
method. The quadratic programming technique (NLPQ), and the cascade strategy converged to design 
with a optimum weight of 951.9 lb which is less than one percent heavier than the optimality criteria 
solution. A lighter design was produced in the proration method that changed both the depths. Design 
became heavier when only the depth of the support member was changed. In other words, intuition was 
not satisfactory for the determinate problem. In summary, the maximum variation in the weight at 7.5 
percent was considered acceptable. Deterministic design methods produced (29 and 33) percent 
variation in the design variables (di and d 2 ), respectively. In other words, the variation in the depths was 
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substantial for the simple problem. For all design points, the displacement was the active constraint, 
while the two stress constraints were passive. 

6.4.2 Stochastic Design 

For stochastic optimization, normal distribution was considered for all the primitive random 
variables, see table 1 . Design solutions were generated in the range of: one failure in two samples (or N 
= 2) up to one failure in two million samples (also in the opposite spectrum not listed). Weight versus 
reliability was plotted in figure 10 to obtain the inverted S- shaped graph. 

For the problem, the weight was 802 lbf for one failure in two samples and it is fifteen percent lighter 
than the optimality criteria solution. The weight was fifty percent higher when reliability was increased 
to one failure in two million samples. As mentioned earlier the weight will depend on the nature of the 
assumptions for the primitive random variables listed in table 1 . There may be inaccuracy in the mean 
value and standard deviation for the variable listed in table 1 . At this early stage of the research, it is 
suffice to say that the stochastic approach produced a design for a specified reliability. The solution 
depended on the assumptions and it should be used with caution. 

6.4.3 Nature of Inverted S-shaped Graph 

The nature of the inverted S-shaped graph was explained considering stress as a random variable with 
a normal distribution as shown in figure 11. For one failure in two samples or reliability p = 0.5, stress 
had a mean value of 25 ksi and for the situation weight was 802 lbf. Consider next a more reliable 
design that correspond to one failure in 100, 000 samples. The value of stress increased to 35 ksi from 
25 ksi. Design became heavy at 1141 lbf because of higher induced stress level. Consider next a very 
reliable design with no-failure. For the situation with reliability approaching unity, the magnitude of 
stress can approach infinity and the design weight would asymptotically reach infinity. In other words, 
weight would increase when risk is reduced. The same logic can be extended to the inverted -S graph in 
the low reliability range (0 - 0.5). When reliability approaches zero, the value of the stress would 
asymptotically decrease to zero with very little weight would be required for the structural design. 


6.4.4 Clamped Beam 

The indeterminate clamped beam shown in figure 5 was made of aluminum with properties 
identical to the tapered beam discussed earlier. It had a length of 4a = 192/h. The uniform thickness of 
the beam was set to (t = 2 in.) . The beam was subjected to three transverse load components as shown in 
the figure with a total magnitude of(P = 100/c/p) . Allowable stress limit was (<r 0 = 20 ksi) and tip 
displacement was not to exceed one quarter of one inch (A < 0.25) . The depths (d v d 2 )were the two 

design variables. For stochastic calculations most random variables listed table 1 was used. Few 
additional variables were required. For displacement limitation a mean value of and a standard 
deviation of (ay = 0.05 m.) were used. The mean values for the three load components (at quarter span, 

at mid span and at three quarter span locations were ( jU P] = 20 kip, /u P2 = 45, ju n = 20 ) , respectively, and 
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corresponding standard deviations were (<t pi =5 kip, = a P2 = a P3 ) . The mean values for the two depths 
were (// P1 = 20 kip, p P2 = 45,// P3 =20) , respectively, and corresponding standard deviations were 
( & p\ ~ 5 kip, — cj p2 = (~y p . ) . 


The stress ratio method produced the design for strength: (ct, = 20.2 in.,d 2 = 1 1 ,5/n.) with a weight of 

608.4 lbf. The design converged in five iterations because of indeterminacy. The maximum 
displacement was ( 8 = 0.38 in.) but the requirement was ( 8 < 0.25) . Solutions were generated by different 

methods for both stress and displacement constraints. The proration technique converged to a weight of 
699.1 lbf. The weight by engineering method was 707.5 lbf. The weight calculated by the optimality 
criteria method was intermediate to the proration and the engineering method. The design generated by 
the optimization methods converged to a structure with four percent higher weight than that was 
calculated by the traditional methods. The variation was about fifteen percent for the design variables. 


The formula for the transverse displacement at the center of the beam can be written as: 


8 =138.24 

max 


/ 2d 2 6 + 21 c/^d 3 + d® ' 


dfd| (d? + d 2 3 ) 

The displacement constraint ( 8 may = 0.25 in.) against the beam weight could be graphed to obtain a figure 
that was similar to figure 9 and it was not repeated. The message however is the same. There are an 
infinite number of designs that satisfy the displacement constraint (g = 8 max - 0.25 = 0) . 


( 9 ) 


For the clamped beam stochastic design solutions were obtained considering multiple random 
variables. Designs were generated for one failure in two samples up to one failure in 1.25 million 
samples. The weight varied from 729 lbf to 3033 lbf, which represented a fourfold increase from the 
mean valued design and a very reliable design with one failure in 1.25 million samples. Likewise, a 
fourfold increase was observed in the two design variables. The displacement limitation was the active 
constraint, see equation 9. The weight and design variables became heavier when risk was reduced. 
Weight versus reliability also traced out an inverted S-shaped graph, which was quite similar to figure 3 
and it was not repeated. 

6.4.5 The Raked Wing Tip Structure 

The raked wingtip structure of the Boeing 767-400 extended range airliner, shown in figure 13 was the 
final example. It was made of metallic and composite materials. The available industrial design was 
considered as the traditional solution. The deterministic optimum design and reliability-based solutions 
were generated for the problem. The component had a wing-box type of construction with face sheets 
and web panels made of composite materials with an aluminum honeycomb core. Its sub-components 
such as the root rib, tip rib, and leading edge parts were made of aluminum and a set of titanium 
members were also used. Its finite element model had about 3000 nodes and 17,000 degrees of freedom 
with 3500 different types of elements. Nominal values were used for Young’s modulus, Poison’s ratio, 
shear modulus, density etc., for aluminum, titanium, composite, and honeycomb materials. The structure 
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was subjected to eight load cases. Strain energy induced in the structure was calculated for each load 
case and it ranged in (16 to 100) units in a normalized scale. Load case 6 contained the peak value for 
strain energy of 100 units. The structure was designed for the critical load case 6 and checked against 
other load cases. The concept to design for a load case with maximum strain energy drastically reduced 
the number of calculations in design optimization, and it was subsequently proved to be correct for the 
problem. Because of proprietary nature of the structure, only a limited amount of information can be 
released. 


For the purpose of design calculation, the wingtip was separated into several substructures made of 
metallic and composite members. The members were grouped to obtain a total of 13 design variables. 
For constraint formulation, the structure was separated into 203 groups of elements to obtain a total of 
203 strain constraints. The composite rod elements were grouped to obtain 16 more strain constraints. 
Three translations and one rotation were constrained at a tip node of the structure for load case 6 as well 
as for load case 7 to obtain eight displacement constraints. The design model has a total of 227 behavior 
constraints. Constraint can be imposed on principal strain, on a failure theory, or on a strain component. 
Design optimization was performed using the testbed CometBoards with MSC/Nastran as the analyzer. 
Only a normalized optimum solution could be given because of the proprietary nature of the data. 

The optimum weight was 16 percent lighter than the traditional design. It should be noted that the 
‘traditional’ design might have included a weight penalty to satisfy manufacturing requirement or other 
constraints. Design changed throughout the structure except for variables 12 and 13, which represented 
the minimum metal thicknesses. Maximum reduction was observed for design variable 8 with a 44 
percent change. The thickness reduced to 0.44 in. if the original value was 1 in. In other words, in the 
original structure the entire region that was associated with design variable 8 represented an 
overdesigned condition. The variable 6 increased to 124 percent from its assumed value of 100 percent, 
or the associated region, was under-designed. The design process redistributed the strain field with many 
active constraints. There was little change in the displacement state or frequency value. 

6.4.5(a) Probabilistic Analysis 

Probabilistic analysis is discussed in brief prior to stochastic design optimization. The analysis 
required a geometrical model, a load model, and a material model. The geometrical model was 
constructed based on the deterministic optimum solution. All 13 design parameters of the model were 
considered as random variables. Their mean values were equal to the deterministic optimum design 
solutions. Their standard deviations were set at 7.5 percent of the mean values. The thicknesses of the 
honeycomb were also considered as random variables. Their mean values were equal to that of the initial 
design with a standard deviation of 10 percent of the mean values. Cross-sectional areas of the bars were 
considered to be deterministic parameters and were equal to the optimum solutions. Load model was 
obtained by considering each component as a random variable with a mean value equal to 66.67 percent 
of the deterministic value, with a 10-percent standard deviation. For example, if a deterministic load 
component is 100 lbf, then its random counterpart would have 66.67 lbf for its mean value with 6.67 lbf 
for the standard deviation. Each load component was changed in a similar manner. Material model 
considered the properties like Young’s modulus, shear modulus, Poisson’s ratio and densities, etc. as 
random variables. The standard deviations for all the material parameters were set at 7.5 percent of their 
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mean values. The mean values were set to 105 percent of their deterministic values for the Young’s 
modulus and shear modulus of fabric as well as the shear modulus for the honeycomb. The mean value 
for Poisson’s ratio was set to 100 percent of its deterministic value. 

6.4.5(b) Probabilistic Solution for Strain and Displacement 

The probabilistic solution or strain and deformation was obtained using the MSC/Nastran code and 
the FPI module of the NESSUS software. Four different types of distributions were considered: normal 
distribution, Weibull function, Lognormal, and Gumbel type 1 distribution. The probabilistic analysis 
produced the mean value of strain (in the element number 1 1531) to be 5050 microstrain by all four 
types of distribution functions. Likewise, the standard deviation remained the same at 16 percent of the 
mean value by all four distribution functions. The mean value of displacement (at node 1 1243 in the z- 
direction) was 9.7 in. with a standard deviation of 10 percent by all four distributions. The observation 
that the mean values and standard deviations did not change for different distribution types was along 
the expected line. 

Probabilistic solutions were obtained for the strain and the displacement for four different types of 
distributions for different failure rates.. For a 50-percent probability of failure, or one failure in two 
samples, the mean value and the calculated strain were equal at 5050 microstrain for normal distribution. 
The same is not true for Weibull, Lognormal, or Gumbel type 1 distributions because of a lack of 
symmetry in such functions. Consider next, 1 failure in 1 million samples. The Weibull prediction was 
conservative. It was about 85 percent of the normal distribution for strain as well as displacement. 
Estimates by Gumbel type 1 distribution was on the higher side. For 1 failure in 1 million samples, it 
predicted about 40 to 45 percent higher strain and displacement than that for the normal distribution 
function. 

6.4.5(c) Stochastic Design Optimization 

Stochastic design calculation used the information given for probabilistic analysis along with 
additional data required to formulate failure modes or the design constraints. It was assumed that mean 
value of allowable strain was 25 percent higher than its deterministic limit, while the standard deviation 
was set at 7.5 percent of the mean value. For example, the limits for strain were set with a mean value of 
s microstrain and a standard deviation of s/10. From deterministic optimization calculations, it was 
observed that only the displacement limitation along the z-direction and rotation had some influence in 
the design while other stiffness constraints were quite passive. Thus, for probabilistic design only two 
stiffness constraints were retained. For the design calculation normal distributions were assumed for all 
random parameters. The optimization calculation required continuous running of the code for more than 
5 days, but the execution was smooth and eventless. The stochastic optimum solution was quite similar 
to that obtained for deterministic optimization. There were numerous active strain, displacement, and 
rotation constraints. The normalized optimum weight was set to 100 units for strength and stiffness 
constraints for 1 failure in 2 million samples. This design exhibited nine active constraints consisting of 
eight strain and one displacement limitation. For an active stochastic constraint the first factor in 
Equation (7) (based on mean values) was equal to the second factor (which was a function of standard 
deviation and the phi critical parameter) in magnitude but with opposite sign. The frequency was 15.94 
Hz for the optimum solution. 
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6.4.5(d) The Inverted-S-Shaped Graph in Semi-Log Scale 

The weight versus probability of success produced an inverted S-shaped graph that was similar to 
figure 10 and it was not repeated. Instead it was graphed in figure 14, in a logarithmic scale, or log(/V), 
along the x-axis. This figure could be approximated into two linear segments. At the intersection point 
(SD) both strength and stiffness constraints became active. From the origin to the point SD, only 
strength constraints were active. From point SD onwards, both strain and stiffness constraints were 
active. The variation of weight was as expected. The weight had a mean value for a 50-percent rate of 
success, or for N= 2. The weight increased when failure rate was reduced. The weight was in agreement 
up to 1 failure in 10,000 samples because only strength constraints were active. The weight increased 
when both strength and stiffness constraints became active. 

The distribution of the strain state in the wingtip was illustrated in figure 15 for three designs 
calculated by the methods. The first design was the traditional design that was obtained by the industry. 
This design might have included a weight penalty from manufacturing consideration. The second 
design was the deterministic solution while the third one represented the stochastic design. The strain 
distribution was uneven for the initial design, and strain exceeded 4000 microstrain at a few localized 
locations. The distribution of the strain state was more even for the deterministic as well as the 
stochastic design solutions. 

Calculations used a Red Hat Linux 2.6.9-67.ELsmp O/S, with x86_64 architecture, 2600 MHz, 4 cpu, 
8 GB of memory, and 32-bit numeric format. One static analysis cycle required about 5 CPU seconds. 
The run time increased to 5 1 seconds for dynamic analysis. Stochastic analysis required about 47 min. 
Deterministic optimization required about 39 min. The CPU time for stochastic optimization was 
enormous at 126 to 128 h of continuous calculations. 


CONCLUSION 

A single winner cannot be identified from the three design methods because acceptable solutions 
were produced by the traditional method, the deterministic optimization and the stochastic design 
method. Designs should be generated following all three methods and compared prior to manufacturing. 
The variation in weight in designs generated by the methods can be considered to be modest. The 
traditional method can generate a design with little less weight than the optimization method. Some 
variation was noticed in designs calculated by the methods. Variation may be due to the indeterminacy 
that can influence stress or strain flow in the members. The traditional design methods should use and 
benefit from the simplified sensitivities of the behavior constraints. Such sensitivity can reduce design 
calculations and may have a potential to unify the traditional and optimization methods. Weight versus 
reliability traced out an inverted-S-shaped graph. The center of the graph corresponded to mean valued 
design with a fifty percent probability of success. A heavy design with weight approaching infinity 
could be produced for a near-zero rate of failure. Weight can be reduced to a small value for the most 
failure-prone design. Reliability can be changed for different components of structure. Probabilistic 
modeling of load and material properties and other design parameters remained a challenge. 
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Figure 3. Comparative evaluations of optimization algorithms 
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Figure 4. Design of a rear divergent flap of a downstream mixing nozzle 
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(b) Rear divergent flap. (e) von Mises stress at optimum. 
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Figure 5. Profile optimization of a cantilever beam for stress and displacement constraints 
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(a) Cylindrical shell. 
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(b) Optimum profile along cylinder length. 

Figure 6. Substructure optimization of a cylindrical shell 
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Figure 8. Tapered beam 
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Figure 9. Design of a tapered beam for strength and stiffness limitations 
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Figure 10. The inverted S-shaped Graph for tapered beam 
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Figure 11. Normal distribution for induced stress 
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Figure 12. Clamped beam 
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Figure 13. Model of raked wingtip of Boeing 767-400 extended range airliner 
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Figure 14. Inverted-S graph in log scale 
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Industrial design: 
Weight = 100 units 



Deterministic design 
optimization: 
Weight =84,5 



Stochastic design 
optimization 
for one failure in 
2 million samples: 
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Figure 15. Optimization process more evenly distributed the strain field in the wingtip 


NASA/CR— 201 3-2 1 7748 


41 



REFERENCES 


1. J.D. Guptill, R.M. Coroneos, S.N. Patnaik and L. Berke, “COMETBOARDS User's Manual,” 
NASA TM-4537, 1996. 

2. S.N. Patnaik, J.D. Guptill, and L. Berke, “Singularity in Structural Optimization,” Int. Jnl. for 
Numerical Methods in Engrg, vol. 36, 931-944, (1993). 

3. S.N. Patnaik, T.M. Lavelle, D.A. Hopkins, and R.M. Coroneos, “Cascade Optimization Strategy 
for Aircraft and Air-Breathing Propulsion System Concepts,” J. Of Aircraft, vol. 34, no. 1, 
136-139,(1997). 

4. RPKJMASTRAN, COSMIC, University of Georgia, Athens, GA. 

5. MSC/NASTRAN Quick Reference Guide, Ver. 68, MacNeal-Schwendler Corporation, 1992. 

6. S. Nakazawa: MHOST Version 4.2. Vol. 1: User's Manual. NASA CR-182235-VOL-1, 1989. 

7. V. Vance and V.A. Tischler, “ANALYZE-Analysis of Aerospace Structures With Membrane 
Elements.” Report AFFDL-TR-78-170, Air Force Flight Dynamic Laboratory, Wright-Patterson, 
AFB, Ohio, (1978). 

8. S.N. Patnaik, R.M. Coroneos, and D.A. Hopkins, “Dynamic Animation of Stress Modes via the 
Integrated Force Method of Structural Analysis,” Int. Jnl. of Numerical Methods in Engineering, 
vol. 40, 2151-2169 (1997). 

9. L.A. McCullers, “Aircraft Configuration Optimization Including Optimized Flight Profiles,” 
edited by J. Sobieski, Symposium on Recent Experiences in Multidisciplinary Analysis and 
Optimization, part 1, NASA CP-2327, 1984. 

10. J.L. Klann and C.A. Snyder, “NEPP Programmers Manual,” NASA TM- 106575, 1994. 

11. W.A. Hafez, “Cometnet-User Manual,” IntelliSys, Beachwood, OH, 1996. 

12. W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, “Numerical Recipes Example Book (C),” 
Cambridge University Press, NY, 1987. 

13. S.N. Patnaik, J.D. Guptill, D.A. Hopkins, and T.M. Lavelle, “Cascade Optimization for Aircraft 
Engines with Regression and Neural Network Analysis- Approximators,” AIAA JP, in press. 

14. S.N. Patnaik, R.M. Coroneos, and D.A. Hopkins, “Substructuring for Structural Optimization in 
a Parallel Processing Environment,” Computer-Aided Civil and Infrastructure Engineering Vol. 

15,209-226 (2000). 


NASA/CR— 201 3-2 1 7748 


42 



15. S.N. Patnaik, A. Gendy, L. Berke, and D.A. Hopkins, “Modified Fully Utilized Design (MFUD) 
Method for Stress and Displacement Constraints,” NASA TM-47431, 1996. 

16. Thacker, B. H., Rhia, D. S., Fitch, S. K., Huyse, L. J., and Pleming, J. B., “Probabilistic Engineering 
Analysis Using the NESSUS Software, Elsevier Structural Safety,” Vol. 28, No. 1-2,(2006), pp 
(83-107). 

17. G.N. Vanderplaats, “DOT Users Manual,” Ver. 2.0. Engineering Design Optimization, Inc., 

Santa Barbara, CA, 1989. 

18. A. Belegundu, L. Berke, and S.N. Patnaik, “An Optimization Algorithm Based on the Methods 
of Feasible Directions,” Structural Optimization J., vol. 9, 83-88, (1995). 

19. G.A. Gabriele and K.M. Ragsdell, “OPT A Nonlinear Programming Code in Fortran 
Implementing the Generalized Reduced Gradient Method User’s Manual,” University of 
Missouri, 1984. 

20. L.S. Lasdon and A.D. Waren, “GRG2 User’s Guide,” University of Texas at Austin, 1986. 

21. H. Miura, and L.A. Schmit Jr., “NEWSUMTl A Fortran Program for Inequality Constraint 
Function Minimization — User’s Guide,” NASA CR-159070, 1979. 

22. J.S. Arora, “IDESIGN User’s Manual,” Version 3.5.2, Optimal Design Laboratory, University of 
Iowa, 1989. 

23. IMSL MATH/LIBRARY FORTRAN Subroutines for Mathematical Applications, vol. 3, Chapt. 8, 
903, (1987). 

24. NAG FORTRAN LIBRARY: MARK 15:E04UCF. NAG Fortran Library Routine Document, 
vol. 4, (1991). 

25. S.N. Patnaik, J.D. Guptill and L. Berke, “Merits and Limitations of Optimality Criteria Methods 
for Structural Optimization,” Int. Jnl. for Numerical Methods in Engineering, vol. 38, 3087- 
3120,(1995). 

26. A. Belegundu and P.L.N. Murthy, “A New Genetic Algorithm for Multiobjective Optimization,” 
AIAA Paper 96-4180, Sixth AIAA, NASA, and ISSMO Symposium on Multidisciplinary 
Analysis and Optimization, Bellevue, WA, 1996. 

27. S.N. Patnaik, R.M. Coroneos, J.D. Guptill and D.A. Hopkins, “Comparative Evaluation of 
Different Optimization Algorithms for Structural Design Applications,” Int. J. Numer. Methods 
Eng., vol. 39, 1761-1774, (1996). 


NASA/CR— 201 3-2 1 7748 


43 



28. K.A. Geiselhart, “A Technique for Integrating Engine Cycle and Aircraft Configuration 
Optimization,” NASA CR-191602, 1994. 


29. K.A. Geiselhart, M.J. Caddy, and S.J. Morris, Jr. “Computer Program for Estimating 
Performance of Airbreathing Aircraft Engines,” NASA TM-4254, 1991. 

30. S.C. Sommer and B.J. Short, “Free-Flight Measurements of Turbulent-Boundary-Layer-Skin 
Friction in the Presence of Severe Aerodynamic Heating at Mach Number from 2.8 to 7.0,” 

NASA TN-3391, 1955. 

31. M.J. Caddy, and S.R. Shapiro, “NEPCOMP — The Navy Engine Performance Computer Program, 
Version I,” NADC-74045-30, Apr. 1975. 

32. R.M. Plencner, and C.A. Snyder, “The Navy/NASA Engine Program (NNEP89): a User’s 
Manual,” NASA TM- 105 186, Aug. 1991. 

33. S.N. Patnaik, R.M. Coroneos, and D.A. Hopkins, “Compatibility Conditions of Structural 
Mechanics,” Int. Jnl. of Numerical Methods in Engineering, vol. 47, 685-704 (2000). 

34. R. H. Gallagher and O.C. Zienkiewicz, eds: “Optimum Structural Design”, John Wiley & Sons, 
London, (1973). 


35. S. N. Patnaik, S. S. Pai and D. A. Hopkins, “Precision of Sensitivity in the Design Optimization of 
Indeterminate structures”, International Journal of Physical Sciences vol. 2 (2007), pp.(0 18-032). 


35. S. N. Patnaik, S.Pai, R.M. Coroneos, “Reliability-based design optimization of airframe 
components”, Proc. IMechE, Part G: J. Aerospace Engineering , ( 2009), 223(G7), pp (1019-1036). 


NASA/CR— 201 3-2 1 7748 


44 



Chapter 2 

Optimization Algorithms and 
Cascade Strategy 




Chapter 2 


Optimization Algorithms and Cascade Strategy 


1.0 INTRODUCTION 

Nonlinear programming algorithm or the optimizer plays an important role in design optimization. An 
optimizer is analogous to equation solver in the finite element analysis. Reliable codes are available for the 
solution of a very large number of linear equations. The same cannot be said for optimizers. Most 
optimizer perform adequately for real life design problem with a limited number of independent variables. 
This is because a design problem is nonlinear in nature and require nonlinear mathematical programming 
algorithms. However, if a problem is linear then linear programming would be applicable and a very large 
system can be solved. 

Over the decades several optimizers with computer codes have been developed. Performance is not 
uniform across different algorithms. We have compared and assessed the performance of different 
optimizers. The reliability and efficiency of a number of state-of-the art optimizers were determined from 
their ability to solve a set of problems, also referred to as the testbed. For small problems, the performance 
of most of the optimizers could be considered adequate. For large problems however, the performance was 
adequate only for a few optimizers. No single optimizer could successfully solve all problems in the testbed. 
It was, further realized that improvement to the algorithms, such as a new search direction or a strategy to 
calculate a step lengths, already available in the state-of-the-art optimizers was not likely to alleviate the 
convergence difficulties. For the solution of difficult problems we have devised an alternate approach 
called the cascade optimization strategy. The strategy uses several available optimizers, one followed by 
another in a specified sequence, to solve a problem. A scheme perturbs design variables in a specified range 
between the optimizers. The cascade strategy has been tested successfully in the design of supersonic and 
subsonic airliner configurations and air-breathing engines for high-speed civil transport applications. 

These problems could not be solved successfully without the cascade strategy. 

Solution to a large problem, that cannot be solved by the cascade algorithm, can be attempted through a 
subproblem strategy. In this technique the original problem is replaced by an equivalent sequence of 
subproblems. A subproblem with a few design variables and a small number of constraints can be solved 
with available optimizers. Solution to the larger problem can be attempted by repeating solutions to the 
sequence of modest optimization subproblems until convergence is achieved for the original problem. This 
strategy is applicable to structures and to multidisciplinary systems. For structures, clustering 
substructures generates the sequence of subproblems. For a multidisciplinary system, individual 
disciplines, accounting for coupling, can be considered as subproblems. A subproblem, if required, can be 
further broken down to accommodate subdisciplines. For the convergence of the strategy, adequate 
coupling is required between the subproblems 

- Optimization algorithms and the cascade strategy are the subject matter of this chapter. The following 
chapter IV is devoted to the subproblem solution strategy. 
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2.0 OPTIMIZATION ALGORITHMS 


A design optimization problem can be cast as a mathematical programming problem of operations research 
as follows: 

Find the n design variables {X } within prescribed upper and lower bounds {X u and X' } to optimize (either 
minimize or maximize) an objective function f (X) subjected to a set of 
m i inequality constraints^, (X) < 0 )andm e equality constra\nts(g i (X) = 0) . 

We provide a few basic definitions of the optimization method. 

Design variables {X} represent parameters that can be changed or adjusted. For example, in the design 
optimization of a truss the areas of the bars become the design variables. Thickness becomes the design 
variable for a plate structure and so forth. Design variables for an airframe can be the wing area, its aspect 
ratio and sweep angle etc. For a jet engine pressure ratio and temperature at specified location can be 
Resign variables. 

Objective function f (X)or the function of merit: Weight is the popular objective function in structural 
design optimization. It can be the weight of a truss or the weight of an airframe wing. The objective 
function for an airframe can be gross takeoff weight, mission fuel, NO x emission, and their combination etc. 
jfengine thrust is the popular merit function in a jet engine optimization. 

A constraint g (X) represents a failure mode. For a structural design problem constraint can be imposed 
on the level of stress or a failure theory for a metallic structure, strain in composites, buckling, natural 
frequency, excessive displacement, dynamic instability and other behavioral characteristics. For an airliner 
constraint can be imposed on landing and takeoff field lengths, range of the airliner and excess fuel 
capacity, just to mention a few. For a jet engine constrains compressor discharge temperature and 
pressure, missed approach gradient thrust etc. Additional constraints, also called side constraints can be 
imposed for geometrical feasibility of a design. For example, in a waffle plate structure the spacing of 
stringers or stiffeners should not exceed the thickness. 

In optimization calculation an equality constraint (g (X) = 0) without any consequence can be replaced by 
two inequality constraints as (g(X)< 0 & g(X) > 0) . At the optimum one of the two inequality constraints 
can be active. When an inequality constraint during optimization calculation becomes zero 
(g(X) < 0 => g(X) = 0), then is it called an active constraint, otherwise it became a passive constraint. 
-Optimization calculations works on active constraints, while the passive constraints are sidestepped. 

Feasible region: The set of design points {X}„ that satisfy all the constraints is referred to as the feasible 
region of the design space. No failure occurs inside the feasible region and any designpoint inside the region 
becomes a feasible design. The optimum solution resides in the feasible region. 
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2.1 Search direction f®) and step length 

The feasible space is searched to locate the design point that optimizes the objective function. The search is 
initiated from an initial design point (X} ,n and it moves along the direction {S} in a step (a) _ y^g search 
process can be expressed by the following equation. 

{xr «Ur» = + a _ {s} - (3.1) 

Equation 3.1 is solved iteratively until the optimum solution is reached. The search direction is normalized. 
There are several techniques to generate a search direction. These are described in standard textbooks (ref. 
1-2) and are not repeated here. Few representative methods included: 

Powell's Method: This method uses Hessian matrix [H]to generate the search direction. The Hessian matrix 
is seldom calculated in closed form. It is updated from an assumed identity matrix with the progress of the 
design iteration. It is set to identity matrix when the Hessian is suspected to have been corrupted. For 
structural design problems the Hessian matrix can become singular because of the dependency of stress 
and displacement constraints (ref. 3). 

Steepest Descent Method: The search direction is set to the gradient of the objective function. Though 
simple, its performance is poor. 

Conjugate Direction Method: This can be considered as a modification to the steepest decent method. The 
direction becomes a combination of the steepest decent direction and a fraction of the previous direction. 

Variable Metric Methods: The search direction is generated from the inverse of the Hessian matrix. It is 
updated from an identity matrix and another symmetrical update matrix is added. 

Newton's Method: This method used the function value, the gradient and the Hessian information. 

Method of feasible direction: A direction is defined as a feasible direction provided a small step along the 
direction retains the feasibility of the design. The direction becomes usable provided the objective function 
reduces. 

Step Length ( a ) . Fibonacci technique and the golden search scheme are two popular methods to calculate 
the step length (a) which is the distance of move along the search direction {S} . The new design after the 
move = X + ay) along the search direction must maintain feasibility. The updated or new 
Resign must have reduced the objective function and it must reside inside the feasible domain. 

2.2 The Kuhn-Tucker Conditions 

A design point (X}''° can be considered as a local optimum of a constrained problem provided the Kuhn- 
Tucker Conditions are satisfied. These necessary conditions are derived from the stationary condition of a 
Lagrange function (z (x, A.)) that for inequality constraints can be defined as: 
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(3.2) 


L(X,A) = f(X) + Y J ^ i 9 i 

/= i 

where, A y is the f Lagrange multiplier corresponding to the constraint g y 

The Kuhn Tucker conditions can be interpreted as the stationary condition of the Lagrange function with 
respect to the design variables and the multipliers {X} and can be stated as: 

m 

V x F (X'-°)+£A y V x g(X'-°) =0 

/=i 

> 0 (the Lagrange multipliers are non-negative) 
g(X'"°) < 0(the design point is feasible) 

V x is the vector of partial derivative with respect to design variables {X} (3.3) 

2.3 Optimization Algorithms in CometBoards 

The CometBoards code evaluated the performance of several optimization algorithms (ref. 4). Because of 
nroprietary limitation, the discussion is limited to six algorithms. 

The sequence of unconstrained minimizations technique (SUMT), as implemented in the code NEWSUMT 
(ref. 5), is available in CometBoards. In NEWSUMT, the classical penalty function has been modified to 
improve efficiency and a modified Newton's approach is used to calculate the direction vector while a 
golden section technique is used to determine the step length. 

Sequential linear programming (SLP), as implemented in reference 6 is available in CometBoards. From the 
original nonlinear problem, a linear programming subproblem is obtained by linearizing a set of critical 
constraints and the objective function around the current design point. The linearization process and linear 
solution sequence is repeated until convergence is achieved. 

The method of feasible directions (FD), as implemented in reference 6, is available in CometBoards. In FD, a 
usable feasible search direction is used. A step length along the search direction is calculated by polynomial 
approximation. 

Three implementations of the sequential quadratic programming technique are available in CometBoards. 
These are: SQP of IDESIGN (ref. 7), DNCONG of IMSL (ref. 8), and NPSOL in NAG (ref. 9). In this technique, 
the original nonlinear problem is solved through a sequence of quadratic subproblems. In SQP of IDESIGN, 
_a Lagrange function is approximated. The step length is obtained by minimizing a composite descent 
function. DNCONG of IMSL uses quasi-Newton updates for the Hessian of the Lagrange function while the 
constraints are linearized. The step length for an augmented Lagrange is calculated using a bisection 
method. NPSOL in NAG also uses an augmented Lagrangian. The search direction is generated through a 
quadratic subproblem while step length is calculated using an augmented Lagrangian, which is designed to 
avoid discontinuities as much as possible. 
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The reduced gradient method, RG, as implemented in the code OPT (ref. 10), has been incorporated into 
CometBoards. This method partitions the design variable into decision and slave variables and a reduced 
gradient is used to generate a search direction. A line search is carried out by bounding the minimum and 
then calculating the minimum within some tolerance. 

The optimality criteria method (OC), available in CometBoards, can be considered as a variant of the 
Lagrange multiplier approach. It was developed for structural design problems. In OC (ref. 11), an iterative 
scheme is followed to update the multipliers and the design variables separately. 

3.0 The Numerical Testbed 

Tho numerical testbed of CometBoards includes over 41 structural design problems, most of which were 
taken from me literature. Minimum weightfwas the obiective and a iinkinq strategy was followed to reduce 
the number of design variables. Stress, displacement, and frequency behavior constraints were imposed. 
'Multiple static load conditions and consistent efemental mass tor dynamic analyses weie also considered. 
The load conditions, mass distributions, and behavior limitations were specified to ensure that several 
types of behavior constraints were active at the optimum. The initial design of unity was considered for all 
problems unless otherwise specified. For each problem a consistent set of upper and lower bounds were 
specified. Typically, default optimization parameters and convergence criteria specified in the individual 
algorithms were used. These parameters, however, were changed when convergence difficulty was 
encountered. A brief description of the examples follows. 

3.1 Examples Pla to Pld: A Three Bar Truss 

The popular 3-bar truss, as shown in figure 1, made of steel was subjected to a single load condition. It had 
three design variables, and six constraints, consisting of three stress, two displacement and one frequency 
limitations. The optimum weight was 92.87 lb and one stress, one displacement, and one frequency 
constraints were active. Five optimizers (SUMT, SQP, IMSL, NPSOL and RG) performed satisfactorily. OC 
was inadequate, yielding a 38.6 percent over-design condition. The problem was solved again for three 
different initial designs (the SUMT optimum design, 150 percent of SUMT optimum, and 50 percent of 
SUMT optimum). Results followed the pattern of the earlier test cases for which the initial design was unity. 
Cases (a, b, c, and d) in table 1 refer to different initial designs used to begin the optimization iterations. 
These were unity for initial design, SUMT solution, 1.5 times the SUMT solution and one half of SUMT 
results. 
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3.2 Example P2: A Tapered 10-Bar Truss 

A tapered 10-bar aluminum truss, shown in figure 2; was subjected to two load conditions. It had 10 design 
variables and a total of 25 behavior constraints, consisting of 20 stress, four displacement, and one 
frequency limitations. The optimum weight was 3326.74 lb with 11 active constraints: 8 stress, two 
displacement, and one frequency limitation. Four optimizers, SUMT, SQP, IMSL, and NPSOL converged for 
this example. Optimizer RG failed, and OC was marginal with a 5.6 percent over-design condition. Cray-YMP 
CPU time varied between 1.28 sec for NPSOL and 1.91 sec for SUMT. 
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Figure 2. A tapered ten-bar truss 
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3.3 Example P3: A Tapered Cantilever Beam 


The cantilever truss of example P2, was modeled next using 8 triangular membrane elements, as shown in 
figure 3. The loads and constraints were identical to example P2. The eight thicknesses of the elements 
were considered the 8 design variables. The problem had a total of 21 constraints, consisting of 16 von 
Mises stress, four displacement and one frequency limitations. The optimum weight was 1440.24 lb with 
six active stress constraints. Five optimizers, SUMT, SQP, 1MSL, NPSOL, and RG performed well while OC 
produced a 2.8-percent over-design condition. Cray-YMP CPU time varied from 1.79 sec for IMSL to 11.22 
sec for RG. 


t y 
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Figure 3. A tapered cantilever beam modeled with eight triangular membrane elements 
3.4 Example P4:A 25-bar truss 

A 25-bar aluminum truss as shown in figure 4, had 8 linked design variables, and was subjected to two 
load conditions. It had a total of 86 behavior constraints, consisting of 50 stress and 36 displacement 
limitations. Four optimizers (SUMT, SQP, IMSL, and NPSOL) converged to an optimum weight of 544.73 lb 
with four active displacement constraints. Optimizers RG and OC failed. Cray-YMP CPU time ranged from 
1.64 sec for SQP algorithm to 6.5 sec for optimizer NPSOL. 
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Figure 4. Twenty-five bar truss 


3.5 Example P5: 165-Ft-Tall Antenna Tower 

A 165-ft-tall steel antenna tower with 252 members, as depicted in figure 5, had six linked design variables 
and was subjected to two load conditions. Its overhead dish antenna was modeled as a lumped mass for 
frequency calculations. It had a total of 529 behavior constraints consisting of 504 stress, 24 displacement, 
and one frequency limitations. Three optimizers, SQP, IMSL, and NPSOL, converged to an optimum solution 
of 5299.84 lb with small deviations. At optimum, six stress, 12 displacement, and one frequency constraints 
were active. Optimizers RG and OC failed while SUMT produced a six percent under-design condition. The 
Cray-YMP CPU time varied between 376.83 sec for SQP and 1893.80 sec for NPSOL 
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Figure 5. One-hundred-sixty-five-ft tall antenna tower 


3.6 Example P6: 60-Bar Trussed Ring 

A 60-bar trussed aluminum ring was subjected to three load conditions and had two lumped masses, as 
shown in figure 6. It had a total of 184 constraints consisting of 180 stress, three displacement, one 
frequency limitations. It has 25 linked design variables. The optimum weight was 414.51 lb, and at 
optimum, 22 stress, one displacement, and one frequency constraints were active. Four optimizers, SUMT, 
SQP, 1MSL, and NPSOL, converged to the optimum solution. Optimizer RG failed, whereas OC produced a 4.1 
percent over-design condition with a 1.1 percent constraint violation. Cray-YMP CPU solution time ranged 
from 36.96 sec for SQP to 144.11 sec for NPSOL. 
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Figure 6. A sixty-bar trussed ring. 


3.7 Example P7: Geodesic Dome 

A geodesic dome, shown in figure 8 with a diameter of 240 in. and height of 30 in., was subjected to a single 
load condition. It was modeled using 156 bars and 96 triangular membrane elements. The bars were made 
of steel with modulus E = 30 000 ksi. Membranes were made of aluminum, with modulus E = 10 000 ksi. 
The bar areas and membrane thicknesses were grouped to obtain eight and four linked design variables, 
respectively. The dome had a total of 254 constraints, consisting of 156 stress for the bars, 96 von Mises 
- stress for the membranes, one displacement, and one frequency limitations. The optimum weight obtained 
was 1022.67 lb with 170 active constraints consisting of 168 stress, one displacement, and one frequency 
limitations. Four optimizers, SUMT, SQP, IMSL, and NPSOL converged with small deviations. Optimizers RG 
and OC failed. The Cray-YMP CPU time varied between 448.32 sec for SUMT to 548.36 sec for NPSOL 


NASA/CR— 201 3-2 1 7748 


56 





Figure 7. A geodesic dome made of bar and membrane elements. 

3.8 Examples P8a to P8d: Intermediate Complexity Wing 

An intermediate complexity wing made of aluminum shown in figure 8, was modeled with a total of 158 
elements consisting of 39 bars, two triangular membranes, 62 quadrilateral membranes, and 55 shear 
panels. The elements were grouped to obtain 57 linked design variables. The wing, which was subjected to 
two load conditions, had a total of 321 behavior constraints, consisting of 316 stress, four displacement, 
and one frequency limitations. The optimum design for this problem was obtained from four different 
initial designs, (1) initial design of unity; (2) initial design equal to the SUMT optimum design; (3) initial 
design equal to 150 percent of the SUMT optimum design; and (4) initial design that was infeasible and at 
50 percent lower than the SUMT optimum design. The optimum weight was 387.76 lb and there were a 
total of 119 active constraints, consisting of 117 stress, one displacement, and one frequency limitations. 
For the initial design of unity, optimizers SUMT, SQP, IMSL, and NPSOL reached the optimum within a 3.7- 
percent error margin. Optimizer RG failed to solve the problem. Optimizer OC also failed to converge to the 
optimum solution but produced a 20.1-percent over-design condition. Cray-YMP CPU time varied between 
1075.21 sec for SQP to 8292.35 sec for NPSOL. 


NASA/CR— 201 3-2 1 7748 


57 



Figure 8. An intermediate complexity wing 



4.0 Performance of the Algorithms 

For discussion, the 41 test cases in the testbed were grouped as small, medium, and large problems. The 
number of linked design variables ranged between three and nineteen for the small problems, also referred 
to as Group 1. The problems are designated as Pla to P5, P7, and P9 to P18. The normalized optimum 
weight for small problems obtained by each optimizer is depicted in figure 9. 
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Figure 9. Performance of different optimizers for small problems. 

For medium or Group II problem, the number of linked design variables ranged between 20 and 39. There 
were 12 medium problems, which were designated as P6, P19 to P23, and P30 to P35. The normalized 
optimum weight for the medium problems obtained by different optimizer was graphed in figure 10. 
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Figure 10 Performance of different optimizers for medium problems. 

Problems with more than 40 independent design variables were referred to as large or Group III problems. 
There were 10 large problems which were designated as P8a to P8d, and P24 to P29. The normalized 
optimum weight for large problems obtained by different optimizer was depicted in figure 11. The 
discussion was separated into five categories: 


1. Convergence to the optimum weight, 

2. Number of active constraints at optimum, 

3. Cray-YMP8E/8128 CPU time required to solve the problem, 

4. Singularity in structural optimization, and 
Default optimization parameters. 
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Figure 11. Performance of different optimizers for large problems 
4 . 1 . Convergence to the Optimum Weight 

The normalized optimum weights for all 41 problems, obtained by the six optimizers are depicted in 
figures 9 to 11 for small, medium, and large problems, respectively. In these figures, unity represents 
optimum weight and more than unity indicates an over-design condition, while less than unity is an 
infeasible solution. For the purpose of comparison, a solution with constraint violation of less than one 
percent and weight which is within one percent of the best feasible design is considered optimum. A design 
is acceptable when the constraint violation is less than one percent and the weight is within five percent of 
the minimum obtained by the eight optimizers. Performance of an algorithm with respect to convergence 
to the optimum solution follows. 

SUMT optimizer converged to the optimum solution for 35 out of the 41examples, that included 17 small, 
nine medium, and nine large problems. The SUMT algorithm failed for four problems. These are: one small 
problem (P5), two medium problems (P21 and P23), and one large problem (P8b). For both medium 
problems, the SUMT solution was more than one percent infeasible. For the large problem, SUMT gave an 
under-designed solution of more than five percent. 


NASA/CR— 201 3-2 1 7748 


61 



Table 1. SUMMARY FOR 41 TESTBED PROBLEMS 


Problem 

Problem description and 

Constraints 

Aefivi 


number 

number of design variables 

specified 

SUMT 

SOP 

IMSL 

NPSOL 

RG 

OC 

Pla 

3-bar mm (3 IDV. ID = 1) 

3S. 2D. IF 

IS. ID. IF 

IS. ID. IF 

IS. ID. IF 

IS. ID. IF 

IS. ID. IF 

‘IF 

Plb 

3 -bar truss (3 IDV, ID -OPT) 

3S, 2D. IF 

IS. ID. IF 

IS. ID. IF 

IS. ID. IF 

IS. ID, IF 

IS. ID. IF 

‘IF 

Pic 

3 -bar truss (3 IDV, ID - I.S x OPT) 

3S, 2D. IF 

IS. ID. IF 

IS. ID. IF 

IS. ID. IF 

IS. ID. IF 

IS. ID. IF 

"IF 

Pld 

3-bar tru« f3 IDV. ID =* 0 } x OP]*) 

3S. 2D. IF 

IS. ID. IF 

IS. ID. IF 

IS. ID. IF 

IS. ID. IF 

IS. ID. IF 

“IF 

P2 

Tapered 10-bar truss (10 IDV) 

20S, 4D, IF 

7S, 2D, IF 

8S, 2D. IF 

8S. 2D. IF 

8S, 2D. IF 

‘IS, IF 

"SS, 2D. IF 

P3 

Tapered cantilever beam (8 IDV) 

16S. 4D, IF 

6S 

6S 

6S 

6S 

6S 

4S 

P4 

25-bar truss (8 LDV) 

5GS, 36D 

4D 

4D 

4D 

4D 

‘2D 

(•) 

P5 

I6S feet tall antenna tower (6 LDV) 

504S, 24 D. IF 

•-7S. I2D, IF 

6S. I2D, IF 

6S. 12D, IF 

6S, 12D, IF 

•if 

(•) 

P6 

60-bar trussed ring (25 LDV) 

180S. 3D, IF 

2 IS. ID. IF 

21S. ID. IF 

21S, ID. IF 

I9S. ID. IF 

* I8S, ID, IF 

■|0S, IF 

P7 

Geodesic dome (12 LDV) 

252S. ID. IF 

I62S. ID. IF 

168S. ID. IF 

150$. ID. IF 

162S. ID. IF 

‘18S.1D 

‘I2S 

P8a 

Intermediate Complexity Wing (57 LDV, ID -1.0) 

316S, 4D, IF 

106S. ID 

IPS. ID. IF 

117S, ID. IF 

75S, ID. IF 

*6S 

‘l9S 

P8b 

Intermediate Complexity Wing (57 LDV, ID =OPT) 

3I6S, 4D, IF 

(») 

106S, ID 

117S, ID, IF 

1 17S, ID. IF 

106S, ID 

*42S 

P8c 

Intermediate Complexity Wing (57 LDV. ID -1.5 x OPT) 

3I6S. 4D, IF 

I09S. ID 

*11 SS. ID. IF 

a 88S. ID. IF 

‘2S 

*7S 

(•) 

m 

Intermediate Complexity Win* C57 LDV. ID ~0 3 x OPT3 

316S.4D. IF 

II7S. ID. IF 

‘99S.1D. IF 

(a) 

* 14$ 

*7S 

fa) 

P9 

10- bar truss (10 IDV) 

20S. 4D. IF 

IS. 2D. IF 

8S. 2D. IF 

8S, 2D. IF 

8S, 2D. IF 

'IF 

‘5S. 2D. IF 

PIO 

10-bar tniM (5 LDV) 

20S. 4D, IF 

6S, 2D. IF 

6S, 2D. IF 

6S, 2D. IF 

6S, 2D. IF 

*3S, IF 

6S. 2D. IF 

Pll 

Stiffened 10-bar truss (18 IDV) 

36S, 4D, IF 

8S 

9S 

8S 

9S 

'IS 

SS 

P12 

Stiffened 10-bar truss (2 LDV) 

36S.4D. IF 

2S 

2S 

2S 

2S 

'if 

2S 

P13 

Cantilever membrane ( 8 IDV) 

16S, 4D, IF 

65 

6S 

6S 

6S 

6S 

4S 

P14 

Cantilever membrane ( 1 LDV) 

16S, 4D, IF 

IS 

IS 

IS 

IS 

IS 

IS 

PIS 

Cantilever membrane 16-quad elements (16 IDV) 

32S.4D. IF 

IS.4D 

(•) 

4D 

*4S 

*5S 

4D 

PI6 

Cantilever membrane 32-quad elements (16 LDV) 

64S. 4D, IF 

IS.4D 

4D 

4D 

‘9S 

*4S 

4D 

PI7 

Cantilever membrane 48 -quad elements (16 LDV) 

96S. 4D, IF 

I3S. 4D 

'13S.4D 

I5S, 4D 

“8S 

‘3S 

I3S.2D 

PI8 


128S. 4D. IF 

32S.4D 

'29S.4D 

3IS. 4D 

*I6S 

fa) 

23S.2D 

P19 

60-bar trussed ring (25 LDV) 

ISOS 

38S 

*35S 

38S 

38S 

'I4S 

40S 

P20 

60- bar trussed ring (25 LDV) 

3D 

ID 

ID 

ID 

ID 

‘ID 

ID 

P2I 

60-bar trussed ring (25 LDV) 

IF 

(») 

IF 

IF 

IF 

‘if 

Tf 

P22 

60-bar trussed ring (25 LDV) 

1 80S. 24D 

28S, ID 

“30S, ID 

27S, ID 

29S, ID 

*20S 

IBS. ID 

P23 


24D.1F 

‘id. if 

ID. IF 

IP. IF 

‘id. IF 

‘IF 

“ID. IF 

P24 

Stiffened 60-bar trussed ring (49 LDV) 

252S 

75S 

75S 

75S 

75S 

75S 

*59S 

P25 

Stiffened 60-bar trussed ring (49 LDV) 

3D 

ID 

*1D 

ID : 

(•) 

ID 

ID 

P26 

Stiffened 60-bar trussed ring (49 LDV) 

IF 

IF 

IF 

IF 

00 

IF 

IF 

P27 

Stiffened 60-bar trussed ring (49 LDV) 

252S.24D 

7SS. ID 

75S, ID 

75S. ID 

75S. ID 

76S, ID 

*17S 

P28 

Stiffened 60-bar trussed ring (49 LDV) 

24D. IF 

ID. IF 

ID. IF 

ID. IF 

‘ID 

ID, IF 

ID. IF 

P29 

Stiffened 60-bar trussed nn£ (49 LDV) 

252S. 3D. IF 

44S. IF 

46S. IF 

46S. IF 

46S. IF 

47S. IF 

*3S 

P30 

Stiffened ring (24 IDV) 

72S 

28S 

28S 

28S 

28S 

27S 

28S 

P31 

Stiffened ring (24 IDV) 

3D 

ID 

ID 

ID 

(») 

ID 

ID 

P32 

Stiffened ring (24 IDV) 

IF 

IF 

IF 

IP 

<»> 

•IF 

IF 

P3J 

Stiffened ring (24 IDV) 

72S, 24D 

2SS, ID 

27S.ID 

28S, ID 

28S. ID 

25S, ID 

“18S. ID 

P34 

Stiffened ring (24 IDV) 

24D. IF 

ID. IF 

ID. IF 

ID. IF 

*2D 

ID. IF 

ID, IF 

P35 

Stiffened rme (24 IDV) 

35a21Mf 

ITS- If 

■12SUE 

,17.s. if 



1ZS» If 

IF 


•Optimum weight obtained differ* by more than 5 percent or constraint violation more than 1 pereent (see ref I). 

IDV: Independent design variable ID: Initial design LDV: Linked design variable OrT: SUMT optimum design D: Displacement constraints F: Frequency constraints 
S: Stress constraints 
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The optimizer SQP of IDESIGN, successfully solved 32 of 41 examples, which included 15 small, 10 medium, 
and seven large problems. This optimizer failed to give a feasible optimum design for three small problems 
(P15, P17, and P18), two medium problems (P19 and P22), and three large problems (P8c, P8d, and P25). 

IMSL optimizer DNCONG successfully solved 37 of 41examples, that contained 17 small, 12 medium, and 
eight large problems. DNCONG optimizer failed to optimize the intermediate complexity wing (problems 
P8c and P8d). 

The optimizer NPSOL successfully solved 25 of 41 examples, which consisted of 13 small, eight medium, 
and four large problems. This optimizer failed (with an infeasible design over one percent) for four small 
problems (P15 to P18); four medium problems (P23, P31, P32, and P34); and four large problems (P8d, 
P25, P26, and P28). It produced more than five percent over-design condition for large problem P8c. 

RG algorithm successfully solved 13 of 41 examples that included seven small, four medium, and two large 
problems. RG optimizer failed for 12 small problems. It also failed for seven medium problems and three 
large problems. The optimizer RG failed with well over 100 percent error in the optimum weight for 15 
problems. 

OC technique successfully solved 16 of 41 examples, which consisted of six small, five medium, and five 
large problems. OC failed for nine small, two medium, and five large problems with an error in the optimum 
weight exceeding five percent, as well as for three medium problems with an infeasible design greater than 
one percent. 

4.2 Number of Active Constraints at Optimum 

The number of active constraints at the optimum for all examples is given in table I. Typically different 
optimizers produced identical numbers of active frequency and active displacement constraints. However, 
the number of active stress constraints generated depended on the optimizer of choice. For example, for 
the geodesic dome problem (P7), the number of active stress constraints produced was 168 by optimizer 
SQP of IDESIGN, 162 by SUMT and NPSOL, and 156 by IMSL. Consider also the set of five examples that 
failed to converge, which produced minimum weights between 3.2 to 12.7 percent over- and under-design 
conditions. These examples produced correct number of displacement and frequency constraints, but failed 
to produce the correct numbers of active stress constraints. The deficiency in the number of active stress 
constraints ranged between three for problem P2 to 42 for problem P8a. For these problems the failure of 
the optimizers could be attributed to their inability to produce the correct number of active stress 
_ constraints. This aspect is also described in the section entitled, "Singularity in Structural Optimization” 
section of this chapter. 
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4.3 CPU Time to Solution 


The normalized CPU times on a Cray-YMP8E/8128 computer were recorded for a set of 14 examples. The 
normalization was with respect to the CPU time required by SQP optimizer, except for problems P8c and 
P8d, which were normalized with respect to SUMT time. The CPU time differed among optimizers as 
shown in table II. Even for a small problem (Pla), normalized CPU time differed from 1.0 for SQP to 22.069 
for RG. For a medium problem (P6), normalized time differed between 1.0 for SQP to 3.899 for NPSOL. For 
a large problem (P8a), normalized CPU time varied from 1.695 for IMSL to 1.116 for SUMT and 1.000 for 
SQP of IDESIGN. We observed that variation in CPU time was small for large problems. 

Table II. OPTIMUM WEIGHT AND CRAY-YMP8e/8128 CPU TIME FOR A SET OF PROBLEMS 
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5.0 Singularity in Structural Optimization 

Singularity in structural optimization has been identified for three situations. 

a. The number of active constraints exceeds the number of design variables. Out of the 41 problems 
the 14 examples listed in table III are prone to this type of singularity. 
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b. Linear functional dependencies among a small number of active stress constraints. This type of 
singularity is suspected to have occurred for some of the examples given in table III. 

c. Linear functional dependencies among a small number of active stress and displacement 
constraints. The identification of this type of singularity by mere inspection may be difficult. 

Table III. Problems with active constraints exceeding the number of design variables 
(Singularity can occur in each of these problems) 
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number 

Description 
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Singularity alleviation can reduce calculations and improve reliability of optimizers. The singularity 
issue is illustrated considering the design optimization of the five-bar truss shown in figure 12. Optimum 
solution is to be obtained for the minimum weight condition for stress and displacement constraints. A 
quadratic programming algorithm was used to solve the problem. The design solutions are graphed in 
figures 13(a) and (b) for two cases: 

Case 1: Singularity is disregarded (see fig. 13(a)). 

Case 2: Singularity is alleviated through a singular value decomposition (see fig. 13(b)). 
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Design case 1, for which independent constraints were not separated, exhibited a wide variation in the 
weight, and it failed to converge. Case 2, which used independent constraints, converged to the optimum 
solution. 



10x10 3 



Number of iterations 


Figure 12. A tapered five bar truss 

Figure 13 (a) Standard Optimization Solution with inherent Singularity 
Figure 13 (b) Optimization Solution When Singularity was alleviated 
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6.0 Default Optimization Parameters 

Default parameters specified in optimization codes, that included, convergence criteria, step length, stop 
criteria, active constraint region, iteration limitations, etc. were used to solve the problems. When an 
optimizer failed to converge, the default parameters were changed in an attempt to successfully solve the 
problem. In the solution of the 41 test bed problems, it was necessary to change the default optimization 
parameters quite often in order to reach the correct solution. On an overall basis, default parameters of 
SUMT, SQP, and IMSL algorithms were adequate for the solution of most problems. Most of the default 
parameters for RG and NPSOL were changed to improve their performances. 

6.1 Remark on Optimization Algorithms 

None of the optimizers available in CometBoards could successfully solve all the problems in the testbed. 
Most optimizers, however, can solve at least one third of the examples. A single winner which can be called 
most reliable and efficient optimizer could not be identified. Overall, three optimizers (IMSL, SUMT, and 
SQP of IDESIGN) scored high marks. For small problems, four optimizers (IMSL, SUMT, SQP of IDESIGN, 
SLP, and NPSOL) satisfactorily solved more than fifty percent of the problems. For medium problems, four 
optimizers (IMSL, SQP ofIDESIGN, SUMT, and NPSOL) produced correct solutions for at least half of the 
problems. For large problems three optimizers (IMSL, SUMT, and SQP ofIDESIGN) were found to be 
reliable and efficient. For large problems, the Cray-YMP CPU time was comparable among the optimizers 
that succeeded. Alleviation of singularity improved the efficiency of an optimizer. 

7.0 CASCADE OPTIMIZATION STRATEGY 

From the comparative evaluation of the optimizers we observed two attributes of the algorithms. None of 
the optimizer when considered individually could successfully solve all the problems in the testbed. Every 
structural problem in the testbed could be solved when attempted by a different optimizer available in the 
testbed. In other words, repeated attempts with different optimizers were found sufficient to solve any 
structural design problems in the testbed. The optimizers were used next to solve two sets of nonstructural 
problems: 

(a) Airliner system optimization 

(b) Variable-cycle multi-mission propulsion engine design 

Even the most robust individual optimizer available in the CometBoards encountered difficulty in 
generating optimum solution to either problem. The difficulty can be attributed to factors such as diverse 
_ design variables and distortion of the constraint space. For example, such problems required combining 
variables with different magnitudes and units of measure, like wing area and engine thrust as well as 
pressure ratios. Diverse constraint types like the takeoff and landing field lengths, compressor 
temperatures, velocities, etc. The complexity is further compounded by the large sequences of optimization 
subproblems that have to be solved to design a variable cycle engine. In brief, constraint formulation and 
design variable formulation available in the CometBoards that successfully alleviated deficiency and 
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worked satisfactorily for structural problems were inadequate to solve airliner and engine design 
problems. Improving the two key ingredients common to most algorithms, such as the search direction and 
the step length, and thereby developing a superior optimizer were seriously considered but dropped. We 
believe that such aspects had been considered by the combined efforts of the developers of the one dozed 
or so optimizers available in the CometBoards testbed. We conceived an alternative approach, called the 
cascade optimization strategy (ref. 12). It derives benefit from the strengths of more than one optimizer. 
The cascade strategy, in other words, uses a number of optimization algorithms, one followed by another 
in a specified sequence. A cascade strategy, for example, can be created by using three optimizers, such as 
SUMT, FD, and NPSOL. In this strategy the problem is solved first by SUMT, and an intermediate optimum 
solution is obtained. The second cycle is begun from the SUMT solution with some perturbation. The 
process is continued with the third optimizer (NPSOL). This cascade strategy was found to be superior to 
using any of the three individual optimizers, especially for multiple-mission engine problems. The cascade 
strategy is illustrated considering two examples: an airliner and a jet engine. 


7.1 Design Optimization of an Airliner 

Design optimization of advanced subsonic and supersonic aircraft concepts have been attempted 
successfully through a "soft coupling” of the Flight Optimization Systems (FLOPS, ref. 13) analyzer to the 
design tool CometBoards. The FLOPS analyzer, through its control and eight discipline modules (weights, 
aerodynamics, engine cycle analysis, propulsion data scaling and interpolation, mission performance, 
takeoff and landing, noise footprint, and cost analysis), can evaluate advanced aircraft concepts and 
formulate their designs as a nonlinear programming problem. Options exist for a number of merit 
functions, such as gross takeoff weight, weight of fuel burned, range, cost, and oxides of nitrogen emissions. 
Free variables for the purpose of optimization include wmg area, wing sweep, wing aspect ratio, wing taper 
ratio, and wing thickness-chord ratio, and thrust or engine size. Important behavior constraints are Mach 
number, altitude, approach velocity, jet velocities, mixed thrust, climb thrust, takeoff and landing field 
lengths, maximum turbine temperature, overall pressure ratio, and bypass ratio for a turbofan. The 
resulting multidisciplinary optimization problem has extremely distorted design space, since both design 
variables and constraints vary over a wide range. For example, an engine thrust design variable (which is 
measured in kilo pounds) is immensely different from the bypass ratio variable (which is a small number). 
Likewise, landing velocity constraint (in knots) and field length limitation (in thousands of feet) differ both 
in magnitude and in units of measure. The difficult nature of the design problem is further compounded by 
statistical, empirical equations and smoothing techniques employed in the FLOPS analyzer. The FLOPS 
analyzer, in other words, can become numerically unstable for some combinations of design variables, 

_ especially for a subsonic airliner. Direct solution of the problem by using the most robust individual 
optimizer available in CometBoards could not provide satisfactory results. However, the application of 
some advanced features and unique strengths of the CometBoards design tool, such as a cascade strategy, 
state-of-the-art optimization algorithms, design variable formulation, constraint formulation, and global 
scaling strategy, successfully solved a number of advanced aircraft design problems. The cascade 
optimization strategy is illustrated here for a subsonic airliner. No single optimizer (e.g., SLP, SQP, FD, 
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SUMT, or IMSL) could provide a satisfactory feasible optimum solution. However, a four-optimizer cascade 
strategy was successful in solving the airliner design optimization problem. The four optimizers used were: 


Sequential Linear Programming: The SLP optimizer, which can provide a quick solution, was used as the 
first candidate in the cascade strategy. The SLP optimizer oscillated rather violently for the first few design 
iterations but produced a converged solution in about 30 iterations, see figure 13. For the problem the SLP 
solution was infeasible and 1380.4 lb heavier than the true optimum takeoff weight. 

Sequential Quadratic Programming: The SQP optimizer was begun from the SLP solution with a 4% 
random perturbation. The algorithm converged to an infeasible solution in about 10 iterations as shown in 
figure 13. This solution was 598.9 lb lighter than the SLP results but heavier than the true optimum by 
781.5 lb. 

Method of Feasible Directions: The FD algorithm was begun from the SQP solution with 1% perturbation. 
The FD optimizer produced a feasible design in about 10 iterations that was it was suboptimal by 738.7 lb. 

Sequential Quadratic Programming: The SQP, which was begun with 1% perturbation from the FD 
optimizer solution, converged in about 25 iterations. It produced a feasible optimum solution of 199 275.6 
lb for the takeoff weight of the subsonic aircraft (which was subsequently verified graphically). 

The four-optimizer cascade strategy successfully solved the subsonic aircraft design problem 


e 


E 

2.5x105 a 



Figure 13. Cascade optimization strategy for design of a subsonic airliner 
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7.2 Wave-Rotor-Topped Engine Design 

Conceptually, the wave rotor replaces the combustor in a conventional air-breathing engine. Wave rotor 
topping can lead to higher engine specific power or more thrust for less fuel consumption. Design 
optimization was carried out for a 47-mission-point (specified through altitudes, Mach number, flow rates, 
etc.), wave-rotor enhanced subsonic gas turbine engine with four ports (combustor exhaust and inlet ports, 
compressor inlet port, and turbine exhaust port). The engine performance analysis and the constraint and 
objective formulations were carried out through a soft coupling of NASA Engine Performance Program 
(NEPP, ref. 14) to the design optimization tool CometBoards. To examine the benefits that accrue from 
wave rotor enhancement, the engine was designed by declaring most of the baseline variables and 
constraints to be passive while considering important parameters directly associated with the wave rotor 
to be active. The active variables considered were rotational speed of the wave rotor, heat added, and fuel 
flow. Important active constraints included limits on maximum speeds on all compressors, 15% surge 
margin for all compressors, and maximum wave rotor temperature. The engine thrust was considered as 
the merit function. The wave rotor engine design became a sequence of 47 optimization subproblems. Only 
the cascade strategy could solve the problem successfully for the entire flight envelope. For the mission 
point (defined by Mach = 0.1 and altitude = 5000 ft), the convergence of the two optimizer (SQP followed 
by FD) cascade strategy is shown in figure 14. The first optimizer (SQP) produced an infeasible design at 
67 060.87-lb thrust in about five design iterations. The second optimizer (FD), begun from the SQP solution 
with a small perturbation, produced a feasible optimum design with an optimum thrust of 66 901.28 lb as 
shown in figure 14. The optimum solution is verified graphically in figure 15. In this figure, observe the 
differences between the individual-optimizer (NEPP) solutions obtained with manual intervention versus 
the cascade results. The cascade solution produced higher thrust than the NEPP. Furthermore, the 
compressor speed was an active constraint in the cascade technique but passive for the NEPP solution. In 
brief, the cascade strategy was successful for the subsonic wave rotor design optimization problem. 
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Figure 14. Convergence of cascade algorithm for a subsonic wave rotor 



Wave rotor heat added, Btu/sec 


Figure 15. Graphical verification for a subsonic wave rotor 
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For the solution of difficult design problems we have devised an alternative approach referred to as the 
cascade optimization strategy. The cascade strategy uses several optimizers, one followed by another in a 
specified sequence, to solve a problem. A pseudorandom scheme perturbs design variables between the 
optimizers. The cascade strategy has been successfully tested out in the system design of supersonic and 
subsonic airliner and air-breathing jet engines for high-speed civil transport application. These problems 
could not be successfully solved by an individual optimizer. The cascade optimization strategy, however, 
generated feasible optimum solutions for both airliner and jet engine problems. 

8.0 CONCLUSIONS 

Reliable optimum solutions for structural problems can be obtained through individual optimizers by using 
constraint and design variable formulations. Individual optimizers, however, were found to be inadequate 
for difficult airliner and engine design problems. A cascade strategy designed by combining a number of 
robust optimizers successfully solved several subsonic and supersonic airliner problems, multimission 
high-speed civil transport engine design problems, and wave-rotor topped engine design problems. The 
open issues with the cascade strategy include the sequencing of optimizers, the individual stop criterion, 
and the pseudorandom perturbation for specific problems. Careful planning and strategizing of the issues 
can provide a successful cascade strategy that will be robust and numerically efficient. 
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Chapter 3 


Subproblem Solution Strategy 


1.0 INTRODUCTION 

Solution of a large design optimization problem requires a robust nonlinear mathematical programming algorithm 
of operations research. Chapter 2 was devoted to different algorithms and the cascade strategy. 

A robust nonlinear mathematical programming algorithm of operations research is required to solve a large 
design optimization problem. Even the cascade strategy can struggle to solve a large problem. 

For solution of such a problem, a subproblem solution strategy has been developed. In this strategy, the 
large original problem is divided into a sequence of modest subproblems that can be solved using available 
algorithms. Solution to the large problem can be obtained iteratively through repeated solutions to the 
subproblems. Subproblem strategy has been implemented in the design testbed CometBoards in sequential as 
well as in parallel computational environments. Several issues were encountered during subproblem solution. 
The issues and their resolutions are discussed in this chapter. 

The basic concept of transforming a large problem, with many design variables and constraints, into a sequence 
N smaller subproblems (/ < / < jv ) is depicted in figure 1 . Each subproblem becomes a modest optimization 
problem (i) with a few design variables and a small number of constraints. The solution to a subproblem can be 
generated using available nonlinear programming optimization algorithms discussed in chapter 2. The solution to 
the larger problem can be attempted by repeating the solution to the set of N subproblems until convergence is 
achieved for the original problem. Solutions to all subproblems can be obtained simultaneously by assigning 
different subproblems to different processors in an advanced computer with a parallel processing features. 
Subproblem optimization, in other words, can become an attractive proposition for the solution of large problems 
provided that the strategy converges and the subproblems can be solved simultaneously with an efficient 
utilization of the available a processors in multiprocessor computer. 

During the subproblem implementation in the testbed CometBoards, several issues were encountered. These 
pertained to: 

1 . Coupling and constraint formulation 

2. Convergence to different local optimal solutions 

3. Amount of computation 

4. Utilization of processors in a parallel computational platform 

To illustrate the issues, a set of problems was selected. Each problem was solved three times: first, 
using regular optimization in sequential computation without a subproblem strategy. Next using a subproblem 
strategy in sequential computation. Finally using a subproblem strategy in a parallel computational mode. The 
three design solutions were compared to illustrate the issues along with their resolutions. 
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2.0 SUBPROBLEM STRATEGY 


There are four principal steps in the subproblem solution strategy. The steps are illustrated considering, the 
design optimization of a cargo-bay support of the International Space Station, shown in figure 2, as an example. 
The support structure is made of four plates, a cluster of plates referred to as a box, and five beams. 

Step 1 : Subproblem model 

A subproblem for a large structural design problem can be created in two stages. First, the original structure is 
divided into several substructures or components. The components are next clustered together to obtain the 
sequence of subproblems. 

Consider the example of the support system, shown in figure 2. It is divided into five substructures. 

The first substructure is the closed box (FGHKIJ) consisting of five plates. Its finite element model is 
obtained by discretizing it by 72 (QUAD4) shell elements. 

The second substructure is a trapezoidal plate ( FHEC ), and its finite element model has 37 shell 
elements. 

The third and fourth substructures are triangular plates (GHE and GHD) with 12 finite elements each. 

The fifth substructure contains the five connecting beams BD, BG, AD, AG, and AF. The beam designs 
included in the fifth substructure are not changed. They become the passive variables and do not directly 
participate in the design optimization calculations. For the purpose of reanalysis, the passive members 
are discretized by four beam elements each, or a total of 20 elements for the fifth passive substructure. 

The finite element analysis model should be validated for the initial design and for the final design. Two 
independent finite analysis codes LEHOST (ref. 1) and MSC/NASTRAN (ref. 2) are used to verify the model of the 
support structure. The model with a total of 133 QUAD4 shell and 20 beam elements is considered adequate for 
design optimization because both analyzers produced an acceptable level of accuracy for stress, displacement, 
and frequency. Two different designs for the support structure are considered for optimization: one with stiffeners 
referred to as a stiffened model, and another without stiffeners, referred to as an unstiffened model. The 
subproblem optimization strategy is illustrated for the unstiffened model. 

The substructures are clustered next to obtain a set of subproblems. A subproblem contains a substructure along 
with some or all of its neighboring substructures. The four substructures are clustered to obtain the following four 
subproblems. 

Subproblem 1: Substructures 3 and 4 
Subproblem 2: Substructures 1 and 4 
Subproblem 3: Substructures 1 and 2 
Subproblem 4: Substructures 2 and 3 
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Notice the coupling between substructures. Substructure 1 is common to subproblems 2 and 3. Likewise, 
substructure 2 is common between subproblems 3 and 4, and so forth. Adequate coupling between substructures 
is required for convergence. Inappropriate coupling can increase the amount of computation and can lead to 
convergence difficulties. Convergence difficulty because of inappropriate coupling is subsequently illustrated in 
this chapter. At this time, substructure coupling is decided intuitively from experience. However, it may be 
possible to automate substructure coupling through the gradient scheme developed in reference 3. 

Step 2: Definition of the original problem 

The design of the large original problem is defined by a set of design variables, an objective function, and a 
number of constraints. For the support system, the thicknesses of different plates are considered the design 
variables. The minimum weight is the objective function. The behavior constraints included Von Mises stress, 
beam column buckling, skin crippling, nodal displacements, and specified frequencies. The design is cast as the 
following mathematical programming problem. 

Find the n-number of design variables {X } within the prescribed upper and lower bounds (x z < X. < X to 
minimize a weight function W f X J under a set of m inequality constraints (gj < 0, j = l,2,...,m) .The weight function can be 
written as 


w=Xp p 

Plates 


f4l> k t k 


V k=l 


+ X Pb L b 

Beams 


( 2 

Z c k b k d k 


V k=l 


\ 

y 


where, 

(p p , p b ) ; weight densities of plate and beam members, respectively 
(t k , A k ) : thickness and associated area of plate 
L b : Length of a beam element 

c k : A coefficient that accounts for the nonprismatic nature of a beam cross section 
(d k , b k ) : depth and width of k tk beam member, respectively 


(4.1) 


The weight along with its gradient are generated in closed form through repeated application of the chain rule of 
differentiation. 

The behavior constraints are specified as follows: 

2.1 Stress constraints 

Constraints to safeguard material yielding or stress limitatiions are specified as 
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g„j = ~i — 1-0 (j = 1 , 2 ,..., ml) (4.2a) 

a jo 

where, is the yield strength 

is von Mises’s yield stress and it is defined as: 

For a shell element: a k = ^a x + a 2 - a x a y + 3(r xy + + t 2 x ) 

For a beam element: a k = ^cr + 3(x 2 y + t 2 x j 

Normal stresses are: ct x and a, and shear stresses are: T xy ,T yz ,T zx 

Stress constraints are imposed on groups of plate or beam elements referred to as a design patch. All stress 
constraints within a design patch are rearranged with the maximum violated constraint at the top and the least 
violated constraint at the bottom of an array. Only a few critical constraints from the top of the array are included 
for the purpose of design optimization. 

22 Buckling constraints 

Constraints that safeguard against buckling under a combined axial force and a bending moment are generated 
by using a beam-column model and are specified as 

Sbj = | F bj | - 1 - 0 (j=m 1 +l,m 1 +2,...,m 2 ) 

where F b represents a critical compressive stress ratio. (4.2b) 

Buckling constraints grouping strategy followed the stress constraint grouping scheme. 

2.3 Crippling constraints 

Constraints to safeguard against crippling are specified as 

§bj = | F cj | “ 1 - 0 (j = m 2 +l,m 2 +2,...,m 3 ) 

where, F c represents a critical cripling ratio. Cripling constraint are also grouped like the stress 

and buckling constraints. (4.2c) 

2.4 Displacement constraints 

Nodal displacement limitations are specified as 

g uj = — --1<0 ( j = m 3 +l,m 3 +2,...,m 4 ) (4.2d) 

U j0 

where, u . , u j0 represent the displacement value and its allowable limit, respectively 
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2.5 Frequency constraints 

Frequency constraints are imposed as 


-^° (j=m 4 +l,m 4 +2,...,m) (4.2e) 

where, f . and f mjo represent the frequency and its allowable limit for the m ,h vibration mode, respectively. 


Step 3: Definition of design optimization for a subproblem 

A subproblem is formulated from the original problem by separating the design variables and constraints 
into active and passive sets. The design variables for a subproblem are specified during its definition and are 
considered active, whereas all other variables of the original problem are considered passive. For example, the 
design variables of the second subproblem associated with substructures 1 and 4 become active design 
variables, and its passive variables include those for substructures 2 and 3. The merit function associated with the 
active design variables become the objective function for subproblem optimization. 

The constraint formulation is carried out in two steps. In the first step, a small set of critical constraints typically 
equal to the number of design variables is selected from among all the specified behavior constraints. The critical 
constraint list should include all types of behavior constraints. In the case of the support design, the constraints 
should include stress, buckling, crippling, displacement, and frequency limitations. In the second step, which is 
specific to subproblem strategy, the constraints are separated into local and global sets. Stress, buckling, and 
crippling within a subproblem are considered local constraints, whereas displacement and frequency are 
considered global constraints because their evaluation requires the response of the entire structure as a single 
unit. The constraint for a specific subproblem is obtained by including all its local and global constraints. 

Constraint formulation also depends on the designer's experience and intuition. Convergence can be achieved for 
appropriate constraint formulation; otherwise, difficulties can be encountered. 

Step 4: Subproblem solution strategy 

The subproblem solution strategy for sequential and parallel calculations is depicted in the flow diagrams of figure 
3a and figure 3b, respectively. The flow diagrams shown in figure 3 have two loops. The inner loop (I = 1 , 2,... N) 
is associated with the optimization of N subproblems. In parallel computational environment, the subproblems are 
assigned to different processors. The outer loop repeats the solution of the N subproblems several times, which 
is referred to as cycles, until convergence is achieved. In the sequential algorithm, the subproblems 
are optimized in sequence, and the design variables for the entire structure (also referred to as the global 
variables) are updated as soon as a solution to any single subproblem is available. For example, the global 
design variables for the entire support system, which is divided into four subproblems, are updated four times in 
four intermediate steps after the solution of each of the four subproblems (see figure 3a). In parallel computation, 
wherein subproblems are assigned to four different processors (see figure 3b), the global design variables can be 
updated only after all solutions have been completed for all four subproblems. In other words, the sequential 
algorithm that benefits from intermediate improvement to the global design variables can converge faster than the 
parallel algorithm because updating the global variables in this scheme cannot proceed until a full subproblem 
optimization cycle has been completed, i.e., after the solution to all four subproblems for the cargo-bay support 
system has been completed. 
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3.0 NUMERICAL EXAMPLES 


The seven problems have been selected to illustrate the difficulties that may be encountered in subproblem 
optimization. The first problem is presented in detail. Results are summarized for other problems. 

Problem 1 : Profile optimization of a beam 
Problem 2: Design of a cylindrical shell 
Problem 3: Design of a cargo-bay support 

Problem 4: Design of spacer structure for International Space Station 
Problem 5: Design of a trussed ring 
Problem 6: Design of a geodesic dome 
Problem 7: Design of a three bar truss 


3.1 Problem 1: Beam profile optimization 

The fixed beam shown in figure 4 is considered to illustrate the difficulties that can be encountered because of 
inappropriate coupling and convergence of substructure strategy to different local solutions. 

The symmetrical steel beam (with Young's modulus E = 30,000 ksi, mass density p = 0.289 lb/in 3 , and length L = 

240 in) is subjected to a concentrated load of magnitude P = 10 kips at its center. The initial depth and width of 
the beam are 12 and 1 in, respectively. The problem is to determine the optimal depth profile for the beam within 
the maximum and minimum bounds of 12 and 1 in, respectively, for the minimum weight condition under a 
combination of stress, displacement, and frequency constraints. The problem uses COSMIC NASTRAN (ref. 4) as 
the analyzer. The finite element model of the problem discretized by 16 beam elements is considered adequate 
because this model produced an acceptable level of accuracy for stress, displacement, and frequency 
constraints. The optimal depth profile for the beam for a uniform width of 1 in is obtained by selecting the depth at 
different locations and linearly interpolating it between the locations. The specified behavior constraints are (1) 
one Von Mises stress constraint (see equation 4.2a) with a permissible strength of 20,000 psi is imposed for each 
beam element, (2) a single displacement constraint is imposed at the mid-span of the beam (see equation 4. 2d) 
with an allowable value of 0.45 in, and (3) a frequency limitation at 40 Hz for the fundamental mode is specified, 
as given by equation (4.2e). Symmetry is exploited, and only half the beam with eight elements is used for design 
optimization. For each element, the depth is considered as the active design variable, and the thickness is 
considered the passive variable. The problem has nine active design variables. The sequential quadratic 
algorithm NLPQL (ref. 5) is used as the primary optimizer, which is supplemented by the optimizer, SUMT or 
Sequential Unconstrained Minimization Techniques (ref. 6) when required. The problem is solved on an RS6000 
IBM workstation. 

The optimal design is obtained for the following seven constraint combinations: case 1: stress; case 2: 
displacement; case 3: frequency; case 4: stress and displacement; case 5: stress and frequency; case 6: 
displacement and frequency; and case 7: stress, displacement, and frequency. For illustration, two arbitrary 
optimization models are considered: The first model has six subproblems, and the second model has three 
subproblems, as follows: 

3.1.1 Optimization Model 1 

Subproblem 1: 4 beam elements (1 through 4) 

Subproblem 2: 5 beam elements (1 through 5) 

Subproblem 3: 5 beam elements (2 through 6) 
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Subproblem 4: 5 beam elements (3 through 7) 

Subproblem 5: 5 beam elements (4 through 8) 

Subproblem 6: 4 beam elements (5 through 8) 

3.1.2 Optimization Model 2 

Subproblem 1: 4 beam elements (3 through 6) 

Subproblem 2: 4 beam elements (1 through 4) 

Subproblem 3: 4 beam elements (5 through 8) 

The second model was used only when the first model experienced convergence difficulty, such as for case 1 
(stress constraints) and case 4 (stress and displacement constraints). For the problem, Von Mises stress is the 
local constraint, whereas both displacement and frequency are the global constraints. For example, subproblem 3 
of model 1 has five local stress constraints, one for each of the five beam elements. Displacement and frequency 
are the two global constraints of the subproblem. 

The results obtained for the beam by using regular optimization and subproblem optimization strategy are 
summarized in Tables 4.1 to Table 4.4 and are illustrated in figures 4, 5 and 6. For the problem, a uniform design 
is also generated by linking all nine design variables to a single uniform-depth variable. The uniform beam also 
was solved for the seven constraint cases, and these results are given in Table 1 . 

From the information given in Tables 4.1 to Tables 4.4 and figure 4, we observe the following. 

The beam profiles obtained for model case 1 (stress) and case 2 (displacement constraints) are similar, as seen 
in figure 5b and figure 5c, respectively. For case 1 , the optimal beam profile represents a fully stressed condition 
because all eight stress constraints are active. The profile has a minimum depth at the quarter span, where the 
bending moment would be zero. The profile peaks toward the center and at the fixed ends, as expected. The 
profile for the displacement constraint is similar to that for the stress constraints except for minor deviations. 
However, the beam profile obtained for the frequency constraint (case 3) is distinctly different from that obtained 
for the stress and displacement constraints (see figure 5d).For the frequency constraint, the optimization process 
attempted to eliminate the mass of the beam at the center-span because this location corresponded to the peak 
amplitude or high accelerations. 

3.1.3 Substructure coupling 

For the beam example, the six subproblem model, or optimization model 1 , was used for the solution of all seven 
cases. However, only five of the seven cases could be solved successfully with this optimization model. For cases 
1 and 4, model 1 was not successful. Case 1 and case 4 were solved successfully by using the second model 
with three subproblems, or half that of model 1 . In the subproblem strategy, an increase in the number of 
subsproblem (as in optimization model 1) need not lead to a successful solution. 

3.1.4 Uniform-depth versus variable-depth beam design 

For all seven cases (cases 1 to 7), the uniform-depth design is heavier than that for the profiled beam, as 
expected, see Table 4.1 . The overweight for the seven cases was between 16 and 90 percent. The maximum 
overweight condition of 90 percent occurred for case 3, or for the frequency constraint. The minimum overweight 
condition occurred for case 2, or for stress and displacement constraints. At optimal design, the uniform beam 
produced one active constraint for five cases and two active constraints for two cases, as shown in Table 1 . The 
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profiled beam, however, produced numerous active constraints. The uniform-beam design required a small 
amount of CPU time that ranged from 7 to 18 percent that required for the variable-depth design (see Table 1). 

3.1.5 Local optimal solutions 

The subproblem strategy can converge to local solutions that can be different from those obtained when regular 
optimization is used, see Table 4.2, figure 5 and figureb 6. However, it was observed that the variation in the 
minimum weight between the solutions was small, ranging from 0.0 to 0.3 percent. To some extent, the mild 
variation in weight can be attributed to numerical accuracy and convergence parameters of the algorithms. In 
other words, the minimum weight remained about the same for the subproblem and regular optimization 
strategies. There was substantial variation in the design variables, or depths of the beam for regular and 
substructure optimizations. The maximum deviation from the regular and subproblem optimizations for beam 
depths for the seven cases was between 10 and 15 percent. Consider case 7, which was the design for stress, 
displacement, and frequency constraints. Beam profiles obtained for subproblem and regular optimizations can be 
separated into three regions: near fixed ends, around the quarter span, and near the mid-span, see figure 5. First, 
the depth peaked near the fixed end with values of 13.4 and 14.6 in for regular and subproblem strategies, 
respectively. Second, the depth variations at the quarter span location were 3.6 and 1 .8 in for regular and 
subproblem strategies, respectively. Third, at the center span, the depths were 10 and 12 in for regular and 
subproblem strategies, respectively. The beam profiles obtained for regular and subproblem optimization 
strategies for four cases (7, 1 , 2, and 3) are depicted in figure 5. The depth variation for all four cases is more 
uniform for regular optimization than for subproblem optimization. For example, for case 7, or optimization with 
stress, displacement, and frequency constraints, at optimum both strategies produced four stress and one 
frequency constraint. Subproblem and regular solutions produced the same set of active constraints. The optimal 
solutions obtained by the regular and subproblem strategies for the other six design cases (1 through 6) followed 
the pattern for case 7 with some variations. The optimal depths for the cylindrical shell (discussed under problem 
2) obtained using regular and subproblem strategies are depicted in figure 9. The profile obtained using regular 
optimization is more uniform than that generated from subproblem optimization. The two optimal weights at 
1161.95 and 1154.1 lb between regular and subproblem strategies, which differed by 0.676 percent, can be 
considered negligible, especially for the problem complexity. The depth differed substantially between the two 
solutions. At the crown, the optimal depths at 1.322 in. and 2.471 in. varied by 53 percent. At optimum, the regular 
optimization and subproblem strategies produced a different number of active constraints. 

3.1.6 Computation time 

The normalized CPU time (on a RS6000 IBM workstation) required for the seven design cases for problem 1 by 
the regular and subproblem strategies is given in Table 4.3. For these problems, the time ratio (the ratio of the 
CPU time for the subproblem optimization scheme to that for the regular optimization) varied from 3.9 for case 2 
(displacement constraints) to 10.8 for case 5 (stress and frequency constraints). The average CPU time ratio for 
all seven design cases was 6.6. On an average for these cases, the subproblem strategy is 6.6 times more 
expensive than regular optimization. A brief description of the other six examples solved using subproblem 
strategy is given next. 

3.2 Problem 2: Cylindrical shell 

The cylindrical shell with rigid diaphragms under two line loads is shown in figure 7a. It has a radius of 100 in and 
a length of 200 in. The minimum weight design for the shell is obtained for stress and displacement constraints. 
Because of symmetry, only one-eighth of the shell is considered for optimization. A finite-element model with 1 00 
shell elements, considered adequate to predict its response, is used for its design. The one-eighth shell is divided 
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into four substructures along its length. The substructures are clustered to obtain three and four subproblems for 
sequential and parallel computations, respectively. An additional fourth subproblem was considered for parallel 
calculations because this scheme experienced convergence difficulty with three subproblems. A single design 
variable is considered for each substructure through profiled depth link factors to provide a heavier design under 
the load. 

The optimal solutions obtained are given in Tables 4.2 to 4.5. The optimal weights obtained by the three 
strategies are comparable at 1 161 .95 lb for the regular optimization, 1154.10 lb for the sequential calculations, 
and 1 153.03 lb for parallel computations. The maximum deviation in the optimal weight for three schemes was 
0.78 percent. However, the design variables differed substantially. The average depths obtained by the three 
schemes at 1.02, 1.15, and 1.15 in, respectively, exhibited a deviation of 11 percent between the regular and 
subproblem strategies. The three maximum depths of 1.32, 2.47, and 2.47 in, respectively, exhibited a variation of 
47 percent. The minimum depths at 0.8344, 0.6879, and 0.6877 in, respectively, differed by 21 percent. The 
depth of the shell under the load along its length is shown in figure 9. Like the beam example, the regular 
optimization produced a more uniform variation in the shell depth, whereas the depth variation was nonuniform for 
the subproblem strategy. The peak variation occurred at the center span. The substructure strategy produced 
only active stress constraints, whereas both stress and displacement constraints were active for the regular 
optimization. To verify the existence of two different designs (one obtained from the subproblem strategy and the 
other from the regular optimization scheme), we resolved the regular optimization case by setting the initial design 
equal to the optimal solution that was generated from the substructure strategy. The solution converged to the 
initial design, illustrating that the problem has two local optimal solutions with about the same minimum weight. 


3.3 Problem 3: Cargo-bay support 

The design of the cargo-bay support shown in Figure 3 and mentioned earlier is used as the next example. The 
support system is made of aluminum (Young's modulus E D 9:9_1 03 ksi, Poisson's ratio _ D 0:303, the weight 
density _ D 0:098 Ib=in3, and the allowable strength _o D 30 ksi). Critical design loads were generated from a 
variety of shuttle accelerations and maneuvers 2 i by using a jiite-element model with 9658 elements and 7439 
nodes. The optimal weight obtained by the three strategies, regular optimization, and substructuring in sequential 
calculations and in parallel mode is approximately 34 lb, but the designs differed. The variation in design variables 
was marginal at 0.3 percent for the parallel and sequential computations, but it peaked at 21 percent between the 
substructuring and regular optimizations. 

3.4 Problem 4: Spacer structure for International Space Station 

The spacer structure for the International Space Station is shown in Figure 6b. The aluminum spacer structure is 
fabricated of tubular members. The design is obtained for pseudostatic load conditions generated from a shuttle 
acceleration of 6.75g along the x axis, 2.25g along the y axis, and 6.75g along the z axis (g is the gravity 
acceleration). Stress, displacement, buckling, and frequency constraints are considered. The structure is modeled 
by using 262 beam elements. For the purpose of design, the structure is divided into four substructures that are 
clustered to obtain four subproblems. For the problem, the optimal solution obtained by the three different 
strategies converged to the same design, as given in Table 2. The optimal weights obtained for regular, 
sequential, and parallel schemes are 307.95, 307.87, and 307.88 lb, respectively. 
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3.5 Problem 5: Trussed ring 

A 60-bar aluminum trussed ring is shown in Figure 6c. The minimum weight for the problem is obtained for stress, 
displacement, and frequency constraints under three load conditions. The problem is solved by using two 
subproblems for parallel and three subproblems for sequential computations. The strategies converged to about 
the same optimal solution of approximately 429 lb (see Table 2). The number of active behavior constraints 
differed: eight, six, and eight for regular, sequential, and parallel calculations, respectively. 

3.6 Problem 6: Geodesic dome 

A geodesic dome (Figure 6d) consisting of 132 bars and 61 nodes is considered as the sixth example. The 
minimum weight design for the problem is obtained for stress, buckling, and displacement constraints. For parallel 
computation, six subproblems were used, whereas for sequential calculations, seven subproblems were used. 
The solution for the problem converged to about the same optimal design with a weight of 92 lb for the three 
strategies, as presented in Table 2. Each strategy produced seven active constraints. 

3.7 Problem 7: Three-bar truss 

A three-bar steel truss, shown in Figure 6e, is designed for the minimum weight under stress, displacement, and 
frequency constraints. The solutions obtained are given in Table 6. This problem is used to illustrate coupling and 
constraint formulation. 

4.0 Issues in Subproblem Optimization 

Issues in substructure optimization are discussed under three topics: (1) coupling and constraint formulation, (2) 
differences in optimal designs, and (3) amount of computation. 

4.1 Coupling and constraint formulation 

Adequate coupling between substructures is essential for the success of the substructure strategy. Inadequate 
coupling can lead to convergence difficulties. Coupling is elaborated further considering the three-bar truss shown 
in Figure 6e. Its three-bar areas are considered as the design variables. The truss minimum weight is the design 
objective. Stress, displacement, and frequency are the behavior constraints. Two different cases referred to as 
uncoupled and coupled subproblems are considered. The uncoupled strategy considers each bar area as a 
subproblem for design, but the entire structure is analyzed for generation of the behavior constraints. In this 
strategy, there are three uncoupled subproblems with no physical overlap between the subproblems except 
connectivity at the free node 4 (Figure 6e). The coupled strategy divides the truss into three subproblems: 
Subproblem 1 contains bars 1 and 2; subproblem 2 contains bars 2 and 3; and subproblem 3 contains bars 3 and 
1 . The subproblems are coupled; for example, bars 2, 3, and 1 are common between subproblems 1 and 2, 2 and 
3, and 3 and 1, respectively, as shown in Figure 6f. For the uncoupled subproblems, a pseudorandom damping 
technique was used to perturb and update the design variables between subproblem solutions. Convergence 
occurred for the coupled strategy even without the use of any damping for all four cases of constraint 
combinations, such as (1) stress, (2) displacement, (3) frequency, and (4) for their combinations, or stress, 
displacement, and frequency constraints (Table 6). The optimal weights obtained by using regular and coupled 
Substructuring for structural optimization in a parallel processing environment 223 subproblems converged to the same 
minimum weight for all four cases. For the uncoupled subproblem strategy, convergence occurred for stress 
constraints only. It is likely that the convergence process was assisted by the compatibility condition that coupled 
the three design variables, or bar areas, such as for the three bars, respectively .24 Convergence dif_culty was 
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encountered for the uncoupled strategy for the other three cases: displacement, frequency, and stress, 
displacement, and frequency constraints. With dif_culty, convergence was achieved when displacements 
were the only constraints and a pseudodamping parameter was increased to 70 percent from its normal 
range of 5 to 10 percent. For the other two cases (frequency and stress, displacement, and frequency 
constraints), the uncoupled strategy failed to converge (Table 6). 

For a successful substructure strategy, behavior constraints should be separated into local and global sets. 
Stresses, Euler buckling, and crippling can be considered as local constraints. Displacement and frequency 
should be treated as global constraints because for their evaluation the entire structure has to be treated as a 
single unit. In substructure strategy, behavior constraints for a substructure should include all its local constraints. 
Ideally, all subproblems should include all global constraints. Convergence was achieved when this constraint 
strategy was followed, whereas difficulty was encountered otherwise. 

4.2 Differences in optimal designs 

As observed for the beam, the optimal designs generated by subproblem and regular optimization schemes can 
be different. The difference is discussed under optimal weight, optimal design, and number of active constraints. 
The deviations in the optimal solution for all the problems are given in Table 4.2. 

Optimal weight: The optimal weight for all examples is presented in the bar chart of figure 6. A minor variation 
in the optimal weight generated by subproblem and regular optimization was noticed for all problems. The 
maximum variation in the optimal weight for all examples was within three-fourths of 1 percent. For 8 out of the 
12 design cases, the difference in optimal weight was less than 0.1 percent. For the other four problems, the 
difference was less than three-fourths of 1 percent. The variation in optimal weight obtained by different strategies 
is very small and can be attributed to convergence parameters and computational inaccuracies. 

Optimal design: Considerable variation in the optimal design can be observed for subproblem and regular 
schemes, as shown in figure 6. These three cases (beam design cases 6 and 7 and the cylindrical shell design) 
exhibited a maximum deviation of about 15 percent in the design variable. A difference of 10 percent or more 
was observed for more than half the examples. The minimum and the average deviations in the bar chart of figure 
6 follow the maximum deviation pattern. In other words, regular and subproblem strategy converged to different 
optimal designs with about the same minimum weight. 

Number of active constraints: The number of active constraints obtained from regular and subproblem 
schemes can be different even though 9 of the 12 problems produced an identical number of active constraints, 
as shown in Table 4.2. For beam case 4, regular optimization produced four active constraints (three stress and 
one displacement), whereas the substructure strategy produced six active constraints (five stress and one 
displacement). The cylindrical shell produced eight and four active constraints for the regular and substructure 
schemes, respectively. The number of active constraints also differed in substructure optimization for sequential 
and parallel calculations. For the 25-design variable ring problem, eight constraints were active for regular 
and sequential optimization schemes, but only six were active for parallel optimization schemes. The optimal 
design of an indeterminate structure need not be unique. However, variations in their minimum weight can be very 
marginal. The subproblem strategy that updates different segments of the structure at different solution stages is 
more prone to converge to different design points with about the same minimum weight. The situation can be 
explained through an examination of the compatibility conditions of indeterminate structures. From the viewpoint 
of practical engineering design, the different solutions can be used for advantage, such as, for example, to 
provide clearance at a certain geometric location without a weight penalty. 
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4.3 Amount of computation 

Substructure optimization is computationally more expensive than the regular scheme. The average substructure 
optimization was 6.62 times more expensive than the regular optimization. As compared with regular optimization, 
the number of reanalyses increased by 6.2 times for substructure optimization, and the CPU time followed the 
pattern as expected (see Table 4.3). Computation time in substructure optimization can be reduced by improving 
reanalysis and design sensitivity calculations through substructure analysis and by using approximate concepts 
such as neural networks and regression analysis. Subproblem convergence histories for the beam and the 
cargo-bay support system problem 1 (case 1) and problem 3 are presented in figure 7, and the number of 
convergence cycles is given in Table 4.4. Typically, convergence is achieved in two or three subproblem cycles, 
except for the beam (cases 1 and 4), which required five cycles. The cargo-bay support system in sequential 
calculations converged in two cycles; however, three cycles were required for parallel computations (figure 7b). 
The cylindrical shell problem required five cycles for parallel computation, but only three cycles were sufficient for 
sequential calculations (see Table 4.4). For other problems, the normalized CPU time for the parallel algorithm 
increased from 1.16 for the ring to 1 .74 for the spacer structure (see Table 4.4). In parallel computation, a greater 
number of cycles may have to be executed than those required for sequential calculations. Table 4.5 and figure 
1 0 give the processor utilization and CPU time fractions on the Cray-YMP supercomputer for the five problems in 
parallel computation. The overall performance for the set of problems can be summarized as follows: 

Processor utilization: For the cargo-bay support problem, a 71 percent speedup was obtained; with two 
processors, it peaked at 97 percent, for the geodesic dome. The processor utilization for the five problems 
ranged between 71 to 97 percent, as shown in figure 10 and Table 5. 2. 

Nonparallel time: The time required for executing the serial portions of the code (primarily for read and 
write operations) represented a very small portion (an average of less than 0.1 percent) of the total CPU time 
(see column three of Table 4.5). 

Overhead time: The time required to assign subproblems to processors and to assign common blocks was 
less than 2 percent of the total CPU time (column four of Table 5). 

Idle time: The time due to load imbalance ranged between 1 .6 to 27.8 percent for various problems. Load 
imbalance in substructure optimization occurred because some of the simpler subproblems were solved 
faster than other complicated subproblems. For the geodesic dome problem, the computational load was 
well balanced among the various processors. For the cargo-bay support system design, the load imbalance 
was 27.8 percent because the problem complexity differed among the subproblems (see Table 4.5). 

Speedup time: The speedup time fraction for computation using multiprocessors for different problems (last 
two columns of Table 4.5) ranged between 71 to 97 percent. When the processors are fully utilized, the percent 
speedup ideally should be 100 percent (column six in Table 4.5). When three processors were used, the 
percentage speedup was between 71 and 95 percent; for two processors, between 83 and 97 percent; and 
for six processors, 87 percent. In brief, the speedup time when several processors were used was good for 
all examples. The last column in Table 4.5 provides the multiprocessor speedup time scaled against that for 
the sequential algorithm. Ideally, these times given in column seven should be equal to the times given in 
column six. A difference between the two columns indicates that a penalty has to be paid for parallel 
computation because of inefficient design updates (see figure 4). 
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5.0 CONCLUDING REMARKS 


The implementation of substructure optimization in CometBoards revealed several issues that are discussed in 
this chapter. We do not claim to have resolved all issues satisfactorily, however. Substructure optimization has 
potential, with merits and limitations. As for merits, a large structure can be optimized using substructure strategy. 
Excessive computation is its primary limitation. This deficiency can be alleviated through the development of 
special reanalysis and sensitivity analysis schemes for both local and global constraints. The process would 
reduce computation. Selection of the size of substructures and associated coupling can improve the performance 
of the optimization strategy. This quantification can be attempted through a gradient strategy that has been 
developed for a related research, see reference 3. 

A successful subproblem strategy depends on subproblem coupling and separating the behavior constraints into 
local and global sets. Substructure optimization schemes can converge to different local designs with different 
numbers of active constraints. The local solutions, however, possessed about the same minimum weight. The 
difference in the optimal design was more pronounced for continuous structures, such as for the beam and the 
cylindrical shell examples, than for trusses. Substructure optimization is computationally more intensive than 
regular optimization. Parallel computation can be more expensive than sequential calculation. In parallel 
computation, the individual processors can be used efficiently. 
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Table 4.1 Optimum Solution for a Fixed Beam 
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Table 4.2 Percentage Deviation in Optimal design between regular and substructure strategies 
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Table 4.3 Comparison of CPU time for sequential subproblem solution versus regular optimization 
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Table 4.4 Comparison of computations for parallel versus sequential algorithms 
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Table 4. 5. Summary of results for parallel calculations 
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Table 4.6 Optimal design of three-bar truss for different percent damping 
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Coupling 



Subproblem (1) 

Figure 1. Subproblem and coupling 
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A* 



Figure. 2. Cargo-bay support system of the International Space Station 
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Cycles 



(a) Sequential calculations. (b) Parallel calculations. 


Figure 3. Subproblem solution strategy 
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(a) Fixed beam with uniform depth. 





(b) Profile for stress constraints. 


(c) Displacement constraint. 


(d) Frequency constraint. 



Figure 4. Optimal profiles for a beam for stress, displacement and frequency constraints 
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Figure 5. Different optimal beam profiles with regular and subproblem strategies 


NASA/CR— 201 3-2 1 7748 


100 



16 


14 h 



Problem 1 : Fixed beam (7 cases) 
2: Cylindrical shell 
3: Cargo-bay support 
4: Spacer structure 
5: Trussed ring 
6: Geodesic dome 
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Figure 6. Percent deviation in optimal design for regular and subproblem solutions 


NAS A/CR— 20 13-217748 


101 
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NLPOL optimizer 



_] 

125 


(a) 


25 50 75 100 

Number of design iterations 

Fixed beam case 1 in sequential calculations. 



(b) Cargo-bay support system in parallel computations. *See figure 4. 


Figure 7. Convergence characteristics of subproblem optimization 
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O Lumped mass (0.68 lbf-sec 2 /in.) 
(e) Three-bar truss. 



(f) Three coupled substructures. 

Figure 8. concluded, (e): 3-bar truss (f): Coupled substructures 
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Figure 9. Optimal profile for a cylindrical shell using regular and subproblem strategies 
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Number of processors 

Figure 10. Percentage speedup: multiprocessors versus single processor 
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Chapter 4 


Analyses and Sensitivity Analyses Tools 


1.0 INTRODUCTION 

The test-bed CometBoards is a multidisciplinary design optimization software. It has been used to design structures, aircraft 
and air-breathing propulsion engines. Each design application used a particular analysis tool. For example, structural design 
calculation used finite element analysis (ref.l), aircraft design required the airliner analysis software FLOPS (ref. 2), while 
the NEPP code (ref. 3) became the analyzer for the design of jet engines. An analyzer for a particular application was 
incorporated into CometBoards through a soft coupling feature that allows quick integration of new or used-supplied 
analyzers without any change to the source code. The major analyses codes that have been soft coupled to CometBoards 
included, MSC/Nastran (ref. 1) for structural design, FLOPS for airframe synthesis and NEPP for jet engine design. A set of 
modest analyses codes are also available in CometBoards. These are called DISP, IFM and its variants. DISP is a stiffness 
method code, which was developed by the Air Force. It is popularly known as ANALYZE/DANLYZE (ref. 4). ANALYZE 
representing analysis for multiple static load conditions and DANLYZE represented dynamic analysis. The integrated force 
method (ref. 5) available in CometBoards (ref. 6) is called IFM for code execution. This method is a fmite-element-based code 
wherein forces are treated as the primary unknowns. IFM also replaces the continuum with a finite element discrete model, 
but IFM explicitly constrains the primary force variables {F} to satisfy both the equilibrium equations and the compatibility 
conditions. 


The output from the analyzer is used to form constraints and the objective function of the optimization problem. For example, 


in structural analysis, the stress output is used to formulate the stress constraints like, 


g = 


- 1 . 0<0 


, where the 


induced stress is <T and strength is <T 0 . The number of stress constraints can become too numerous because of the use of 

fine finite element model with many thousand degrees of freedom. There is a need to reduce the number of stress constraints. 
Design optimization becomes computation intensive because of analyses and sensitivity analyses require many calculations. 
The constraint formulation feature of the CometBoards code reduced the number of constraints without any detrimental 
effect. It was noted, that in typical structural design applications, most of the computation, often exceeding ninety nine 
percent of total calculation, was because of the analysis tool (ref. 6). 


Solution to an optimization problem (defined by an objective function and a set of constraints) is generated iteratively in two 
principal steps. At an iteration a search direction (?) is calculated first. Then the design is updated by calculating a step 

length is (X as: ( X new — X old + a ?) , where X new is the new design, X old is the earlier design. The iteration process is 
continued until convergence. 

The calculation of the search direction required the gradients of the constraints and the objective function that is also referred 
to as sensitivity analysis. Several issues pertaining to analyses and sensitivity analyses were encountered during solution to 
multidisciplinary design problems even when well-tested analyses tools were used. Three issues and their resolutions are 
discussed in this chapter. 
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1. Analysis Instability: An optimization algorithm requires values of the constraints and the objective function at any 
arbitrary location in the design space. This information is required throughout the iterative process. If such values 
are not available or a discontinuity is encountered then algorithm will abort abruptly. Design calculations have to be 
reinitiated or the problem cannot be solved. Analysis instability and resolution will be illustrated through solution of 
a real world airliner problem. 

2. Eccentric Design: An optimization algorithm can converge to an irregular local solution that cannot be 
manufactured. Implicit constraints can be the probable cause for the convergence to an eccentric design. An 
augmentation of animation into optimization alleviated the limitation for a specific situation. The issue will be 
illustrated and a technique will be suggested to alleviating the deficiency. 

3. Singularity in Optimization: Structural optimization can become a singular problem because the Hessian matrix 
required to generate the search direction can become rank deficient. The singularity issue and its possible resolution 
are discussed. 

2.0 ANALYSIS INSTABILITY 

Stress and displacement response of very large structures with many thousands of degrees of freedom can be obtained using a 
finite element code like MCS/Nastran. The structural analysis code is robust and it seldom exhibited an instability. This code 
has been soft coupled to CometBoards and a number of problems have been solved. The same cannot be said about a non 
structural analysis code, like FLOPS. The Flight Optimization System (FLOPS) of NASA Langley Research Center is a 
standard aircraft analyzer. The FLOPS code combines multiple disciplines from aerodynamics and engine cycle analysis to 
mission performance. The code uses data tables for internal calculations. For a subsonic aircraft problem the code became 
unstable at several design points. The analysis limitation propagated into design optimization calculations, and the 
optimizer encountered convergence difficulty. The anomalous design points resided in the vicinity of the optimum solution. 
These design points could not be segregated prior to the optimization calculations. The aircraft problem appeared 
to be a good candidate for the application of approximation techniques. To overcome the deficiency two competing 
approximation techniques: neural network (NN) and regression methods were investigated. The regression method used 
a set of basis functions and provided both function and gradient information. NN approximation also used a variety of 
kernel functions and produced the same two pieces of information. Both methods had been applied successfully for a variety 
of multidisciplinary applications (ref. 7-9). The approximate methods were developed using a set of high-fidelity 
training pairs and selected basis functions. The approximate models were validated for use as an alternate reanalysis tool 
for the subsonic aircraft analysis and design optimization. Design optimization of the subsonic aircraft was obtained via 
the CometBoards test bed. The regression method and neural-network based aircraft analysis tools were incorporated into 
CometBoards. The optimum solution of the subsonic aircraft could be obtained using any one of the three analysis methods: 
the FLOPS code NN and regression method based analyzers The design capability was also used to calculate 
sensitivity with respect to the bounds on aircraft constraints such as for example, the takeoff and landing field lengths. 

The discussion is limited to the performance of the three different analysis methods in design of a subsonic aircraft. 

Optimal solutions calculated by the methods were compared. The efficiency in analysis and design was examined 
by comparing the central processing unit (CPU) time to solution. 
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2.1 Subsonic aircraft design problem 


The multidisciplinary FLOPS code can be used for preliminary design evaluation of aircraft concepts. The FLOPS 
FORTRAN code has nine modules: weights, aerodynamics, engines cycle analysis, propulsion data scaling and interpolation, 
mission performance, takeoff and landing, noise footprint, cost analysis, and program control. The FLOPS manual (ref. 2) 
specifies preparation of input data, which follows a name list format with default values. The subsonic aircraft has a fuselage 
with a length of 152.35 ft, width of 16.44 ft, and depth of 17.00 ft. 

The subsonic aircraft was to carry 200 passengers and an eight-member crew, to fly at a cruise speed of 0.8 Mach over a 
range of 2500 n mi. It was powered by two high-bypass-ratio engines with a nominal thrust of about 48 925 lbf. 

It is a separate-flow turbofan engine with two compressor components. The weight of the engine is 9410 lbf. The baseline 
engine nacelle is 19.75 ft long with an average diameter of 7.81 ft. Wing area is 2272 square feet, sweep angle is 31.5°, taper 
ratio is 0.267, and wing thickness-to-chord ratio is 0.109. Nominal parameters of the engine include a bypass ratio of 5, 
overall pressure ratio of 29.5, fan pressure ratio of 1.67, compressor discharge temperature of 1460 °R, and maximum 
dynamic pressure of 800 lbf/sq.ft. Fuel capacity is 57 000 lbf, and there are 10 tanks. The range of the aircraft is 2500 
nautical miles, the maximum cruise altitude is 4000 ft, and the maximum operating Mach number is 0.843. The ramp weight 
is 250 000 lbf. Maximum allowed takeoff and landing field length is 6000 ft. Maximum allowed approach velocity is 125 
nautical miles. Ground operations include a takeoff time of 0.4 min, taxi in and out time of 10 min, and reserve holding time 
of 0.5 hr. 

The objective of the optimization study was to determine the airframe-engine design combination that would meet specified 
constraints and minimize the gross takeoff weight of the airliner. A good match between airframe and engine could be 
achieved by combining the airframe variables with the engine parameters. Nine active variables, listed in Table 1, were 
selected. There were four airframe design variables: wing aspect ratio DV 1? wing area DV 3 , sweep angle DV 4 , and thickness 
to chord ratio DV 5 . The five engine design parameters were: thrust DV 2 , the turbine inlet temperature DV 6 , the overall 
pressure ratio DV 7 , the bypass ratio DV 8 , and the fan pressure ratio DV 9 . Four other variables could be included if required. 
These were: taper ratio of wing DV a , cruise Mach number DV b , cruise altitude DV C and engine throttle ratio DV d . 

The constraint set included: Landing velocity gi not to exceed 125 knots. Field lengths for takeoff g 2 and landing g 3 were not 
to exceed 6000 ft. Missed approach gradient thrust g 4 and second segment climb thrust g 5 were required to be positive. 
Compressor discharge temperature g 6 should not exceed 1460 °R. Excess fuel g 7 should be positive. Constraints g 1? g 2 , g 3 , and 
g 6 restrict the landing approach velocity, takeoff field length, landing field length, and compressor discharge pressure, 
respectively, to not exceed their upper bounds. The g 4 , g 5 , and g 7 constraints, scaled with respect to 101 000, 100 000, and 5 
000 lbf, respectively, restrict the variables to be positive. These are referred to as the positivity constraints. Four additional 
constraints could be included. These were: range of aircraft, specific thrust, specific fuel consumption and compressor 
discharge temperature. 


Table 1 . Design variables and constraints of the subsonic aircraft 


Design Variables (DV) 

Constraints (g) 

Wing aspect ratio, (EAR) 

DVi 

Landing approach velocity, (VAPP) 

gi 

Engine thrust, (ETHRUST) 

dv 2 

Takeoff field length, ( FAROF) 

g2 

Wing area, (ESW) 

dv 3 

Landing field length, (FARLD) 

g3 

Quarter chord sweep angle, (ESWEEP) 

dv 4 

Missed approach gradient thrust, AMFOR) 

g4 

Thickness to chord ratio, (ETCA) 

dv 5 

Second segment climb thrust, (SSFOR) 

g5 

Turbine inlet temperature, (EETIT) 

dv 6 

Compressor discharge temperature, (CDT) 

g6 

Overall pressure ratio, (EEOPR) 

dv 7 

Excess fuel capacity, (EXFUE) 

g7 

Bypass ratio, (EEBPR) 

dv 8 



Fan pressure ratio, (EEFPR) 

dv 9 



Variables considered passive but can be included 

Passive constraints that can be included 

Taper ratio of wing, (ETR) 

DVa 

Range of aircraft, (RANGE) 

ga 

Cruise Mach number, (EVCMN) 

DVb 

Specific thrust, (ST) 

gb 

Cruise altitude, (ECH) 

DVc 

Specific fuel consumption, (SFC) 

gc 

Engine throttle ratio, (EETTR) 

DVd 

Compressor discharge pressure, (CDP) 

gd 
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( 1 ) 


The FLOPS code has a provision to use a composite merit function that can be expressed as 

Obj = L w kPk 

k = 1 

Here, Obj represents the merit function, w k represents the k th weighing factor, and the parameter j3 k can be selected from the 
following list: 

(1) Gross takeoff weight of the aircraft 

(2) Mission fuel 

(3) The product of the Mach number and the ratio of lift-to-drag 

(4) Range 

(5) Cost 

(6) Specific fuel consumption 

(7) NOx emissions 

For the subsonic aircraft problem, the gross takeoff weight was selected as the merit function, by setting w\ = 1.0, and setting 
other weighing factors to zero. The objective of the optimization study was to determine the optimum gross takeoff weight of 
the aircraft for the nine design variables and the seven behavior constraints listed in Table 1. Optimum solution was also 
calculated for the aircraft to operate on shorter and longer runways in the 4500 to 7500 ft range. This study was referred to as 
the sensitivity analysis. The FLOPS code calculated the performance parameters for subsonic aircraft, which in turn was used 
to generate the constraints and merit function required for its design optimization. The code synthesizes eight disciplines: 
weight estimation, aerodynamic analysis, engine cycle analysis, propulsion data interpolation, mission performance, airfield 
length requirements for takeoff and landing, noise footprint calculations, and cost estimation. The FORTRAN code has 1 1 
modules with over 42 000 FORTRAN statements. The subsonic aircraft problem required several input/output (I/O) files. A 
brief description of the code is given in Ref. 2. 

2.2 Justification for Use of Approximate Methods 

Numerical data tables, also referred to as table lookups, used in the FLOPS code can interrupt calculations and the code can 
stop abruptly. Approximate methods can alleviate such limitation of the code. The difficulty encountered in the FLOPS code 
was illustrated by generating the aircraft response for a set of design points that lied in the vicinity of the optimum solution. 
The FLOPS code was run for three sets of analysis data that were created by a pseudo-random perturbation about a base line 
design within prescribed upper and lower bounds as shown in Table 2. The design space spread was modest. It was about 10 
percent on either side of the baseline design. The first set of data was referred to as “small-model,” and it contained 1200 
design points. The “standard-model” and the “large-model” contained 2400 and 4800 points, respectively. Each set of the 
nine design variables and the seven response variables (associated with design constraints) constituted one I/O pair (which 
was also used to train the approximate methods). 

The success rate of the FLOPS analyzer is given in Table 3. The rate of success was about 80 percent for each model. For the 
standard model only 1943 usable I/O pairs could be generated out of the 2400 requested design points. The aircraft weight 
saturated at a quarter million pound- force for 448 design points. The code aborted for seven designs. Turbine entry 
temperature reached a million degrees for one case and a zero thrust condition was encountered for another case. The 250 


NAS A/CR— 20 13-217748 


114 



000 lbf weight and zero thrust condition were either referenced or flagged value of the unsized aircraft. The response for the 
small and large models was similar with minor deviations. 

The design space of an aircraft optimization problem was distorted because both design variables and constraints varied over 
a wide range. For example, an engine thrust design variable measured in kilo pound- force was immensely different than the 
bypass ratio, which was a small dimensionless number. Likewise, a landing velocity constraint in knots and a field length 
limitation in thousands of feet differed both in magnitude and in units of measure. In the design optimization test bed 
CometBoards the effect of distortion was reduced by scaling the merit function, design variables, and constraints such that 
their normalized magnitudes were around unity. 


Table 2. Baseline design variables with bounds 


Design variables 

Lower bound 

Initial design 

Upper bound 

DVi 

Wing aspect ratio 

7.340 

8.500 

8.810 

dv 2 

Engine thrust in lbf 

28000 

31500 

34200 

dv 3 

Wing area in sq. ft 

1830 

2000 

2200 

dv 4 

Quarter chord sweep angle in deg. 

16.0 

18.5 

21 

dv 5 

Thickness to chord ratio 

0.088 

0.095 

0.0997 

$ 

dv 6 

Turbine inlet temperature in °R 

2950 

3000 

3100 

’dv 7 

Overall pressure ratio 

38 

40 

40.50 

sjs 

DV S 

Bypass ratio 

5 

6 

6.10 

*dv 9 

Fan pressure ratio 

1.8 

1.85 

2 


Redundancy in these design variables may cause instability in the subsonic aircraft calculations 


Design optimization of the subsonic aircraft was attempted using the combined CometBoards and FLOPS codes. None of the 
one dozen individual optimization algorithms available in the CometBoards test bed could successfully solve the problem. A 
better solution could be obtained when a cascade strategy was employed. The generation of an optimum solution required 
manual intervention, restarts, as well as a change of the initial designs and bounds. A four-optimizer cascade was employed 
to solve the problem. The optimizers were: sequential linear programming (SLP), followed by a nonlinear quadratic 
programming algorithm (NLPQ), then method of feasible directions (FD), and finally NLPQ. 


NAS A/CR— 20 13-217748 


115 




Table 3. Success rate of the FLOPS analyzer for the subsonic aircraft 


Input/output (I/O) pairs 

Small Model 

Standard Model 

Large Model 

Total number of I/O pairs 

1200 

2400 

4800 

Usable I/O pairs 

991 

1943 

3880 

Bad I/O pairs 

209 

457 

920 

Success rate in percentage 

83 

81 

81 

Saturated at 250 kip for weight 

204 

448 

891 

Code aborted 

3 

7 

16 

Negative one million for engine thrust 

1 

1 

6 

Zero thrust condition 

1 

1 

7 

Used for training 

900 

1800 

3600 

Used for validation 

91 

143 

280 


Solutions generated on an IBM and a SGI workstations were depicted in figure 1 . In the IBM workstation, the cascade 
algorithm converged to a feasible 199 276 lbf for the aircraft weight. The SGI solution (depicted in the lower portion of 
figure 1) converged to a slightly infeasible solution. The same cascade algorithm encountered difficulty on an SGI 
workstation when it was initiated from a different initial design. 
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(a) IBM workstation 



Figure 1 . Convergence history for the subsonic aircraft with FLOPS analyzer and a cascade 
strategy in an IBM and SGI workstations. 


NAS A/CR— 20 13-217748 


117 


Likewise a slightly different cascade exhibited a contrary move resulting in an infeasible solution. The problem was solved 
on the SGI workstation with the original cascade algorithm when the design bounds were changed; see figure 1 (lower 
portion). The optimum designs are given in Table 4. some deviation is observed in the two solutions. There was only a 0.1- 
percent change in the aircraft weight. There was a 0.7-percent deviation in the engine bypass ratio design variable and 3 
percent variation in the second segment climb thrust constraint. Such deviation is considered minor because the subsonic 
airframe engine synthesis is a difficult nonlinear multidisciplinary analysis as well as design problem. The subsonic aircraft 
problem appears to be a candidate for the use of approximation techniques because the FLOPS analyzer can fail at multiple 
design points. The subsonic aircraft optimization process can become tedious. A significant reduction can be achieved in the 
CPU time to solution. 


Table 4. Summary of optimum design solutions for the subsonic aircraft 


Parameters 

IBM Solution 

SGI Solution 

Percent deviation 

Aircraft weight in lbf 

199275.578 

199254.844 

0.10 

Wing aspect ratio (DV1) 

8.547 

8.63 

-0.96 

Engine thrust (DV2), 

31,589.572 

31,595.923 

-0.02 

Wing area (DV3) in sq. ft. 

1897.735 

1879.461 

0.972 

Quarter chord sweep angle (DV4), 

15.650 

16.411 

-4.637 

Thickness to chord ratio (DV5) 

0.093 

0.093 

0.0 

Turbine inlet temperature (DV6) in R 

3060 

3100 

-1.290 

Overall pressure ratio (DV7) 

40 

40.188 

-0.469 

Bypass ratio (DV8) 

5.936 

5.896 

0.678 

Fan pressure ratio (DV9) 

1.824 

1.80 

1.333 

Landing approach velocity (gl) in kn 

119.25 

119.72 

-0.392 

Takeoff field length (g2) in ft 

6000 

6042.66 

-0.706 

Landing field length (g3) 

5490 

5514.84 

-0.450 

Missed approach gradient thrust (g4) 

3737 

3905 

-4.318 

Second segment climb thrust (g5) 

8300 

8548 

-2.901 

Compressor discharge temperature (g6) 

1423.50 

1429.81 

-0.441 

Excess fuel capacity (g7) 

0.2 

0.0 

- 
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2.3 Regression method approximation 

The linear regression method and NN technique were used as two competing approximators in the optimization test bed 
CometBoards. The regression method used several types of basis functions. These functions can be selected from (1) a full 
cubic polynomial, (2) a quadratic polynomial, (3) a linear polynomial in reciprocal variables, (4) a quadratic polynomial in 
reciprocal variables, and (5) combinations thereof. Consider, for example, regression analysis of an ^-variable model 
with a combination of a cubic polynomial in design variables and a quadratic polynomial in reciprocal design variables. The 
regression function has the following explicit form: 


n n n 


y(x) = A,+Z Pi x i +ZZ Pij x i x j + Z X Z Pi 

i= 1 i= 1 7=1 i—\ 7=1 k=\ 


iJkXiXjX k 


n 1 n n 1 


( 2 ) 


/=1 


JC. 


i = 1 7=1 


x,.x. 


The regression coefficients yj3 & j3j were determined by using the double precision general matrix linear least squares 

solver (DGELS) routine of the Lapack library (ref. 10). The gradient or sensitivity matrix of the regression function with 
respect to the design variables is obtained in closed form. For the example with n variables, the gradient matrix for the 
regression function given by equation 2, has the following form: 


Vy = 


d 

cbtj 

d 




d 


dx„ 


( 3 ) 
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where, 


Qy n n-l n n O 

^ = Pt + X /G/ + A< x , + X X Pyf x i x j + X Pm x i x t + A.f x ? — f 

1=1 z=l y=z+l z=l 

P v =Pji for i>j; Pij k =PjJo rj >k , etc. 


J_V R * Pu 

x] tt u x ; . x] 

( 4 ) 


Once the regression coefficients have been obtained from a single training cycle, the reanalysis and 
sensitivity calculations given by equations (2) to (4) required trivial computation. 


2.4 Neural network technique 

The NN approximator Cometnet (ref. 11) is a general-purpose object-oriented library. Cometnet is soft-coupled to the 
CometBoards test bed. The NN capability provides both the function value and its gradient. Cometnet approximates 
the function and its gradient with R kernel functions as follows: 


R n r 
r = 1 z=l 



i’ 5 ") 

( 56 ) 


where y is the functional approximation, x is the vector of independent variables, 
cp ri represent R kernel functions, n r represents the number of basis functions in a 
given kernel, and w ri are the weight factors. 

Cometnet permits approximations by using different types of kernels, which include linear, reciprocal, and polynomial, as 
well as Cauchy and Gaussian radial functions. A Singular Value Decomposition algorithm (ref. 12) for computing the weight 
factors in the approximating function is used to train the network. A clustering algorithm is used to select suitable parameters 
for defining the radial functions. The clustering algorithm, in conjunction with an optimizer, seeks optimal values for the 

parameters over a range for the threshold parameter T within its domain (0 < T < l) . The mean-square error during 

training is reduced by increasing the threshold, which corresponds to an increase in the number of basis functions. Over- 
fitting is avoided with a competing complexity-based regularization algorithm. Different basis functions can be used to 
training the merit function and each of the constraint functions. 
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2.5 Training approximate analyzers 


The I/O pairs generated earlier (see Table 3) were used to train three models for the NN and regression methods. 

The models were referred to as small, standard, and large models. The number of I/O pairs used to train and validate the 
models was (900, 91) for the small model, (1800, 143) for the standard model, and (3600, 280) for the large model. 

Each method had nine free variables, which represented the design variables given in Table 1. Aircraft weight and the seven 
constraints were approximated individually. The basis functions for both approximators contained a full quadratic 
polynomial in the design variables (DV) along with a linear reciprocal expression in the DV. Each approximator had 
64 unknown coefficients. The redundancy (defined as the ratio of I/O pair to number of coefficients) was 14, 28, and 56 for 
the small, standard, and large models, respectively. The values of the coefficients in NN and regression need not be the 
same because they were generated following different procedures. The CPU time for training, reanalysis, and design 
optimization on an SGI octane workstation with the irix 6.5.19m operating system and a 300 MHZ processor was 
given in Table 5. The regression method required a fraction of a CPU second for training. The NN training required 
between 1 and 9 minutes. For a single analysis cycle, the FLOPS code required about 3 CPU seconds. This was 
reduced to milliseconds by the approximators. Gradient calculation was inexpensive by the approximators. For 
optimization the CPU time to solution by FLOPS was 34 minutes. This was reduced to less than two seconds by the 
regression method, while the NN average time was about 4 minutes. For analysis and design calculations the 
approximate methods were found to be efficient. 


Table 5. CPU time in seconds in an SGI octane workstation 


Activity 

Rehression method 

Neural network technique 

Small model 

Standard 

Large 

Small 

Standard 

Large 

Traning, in seconds 

0.2 

0.4 

0.8 

59.1 

136 

538.8 

Re-analysis, in s 
(FLOPS = 3.1 s) 



0.08 



2.4 

Re-analysis with closed 
form gradients, ms 



0.14 



13.5 

Design optimization, s 
(percent of FLOPS 
solution time =2031 s) 

1.6 

(0.78) 

1.7 

(0.84) 

1.6 

(0.78) 

300.9 

(15) 

199.2 

(9.8) 

166.7 

(8.2) 
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2.6 Performance of approximators for analysis and design optimization 


Solutions obtained by different approximation models for a randomly selected design point (DVi, . . . , 

DV 9 = 8.9579, 31607.7515, 2177.9724, 18.5423, 0.0874, 2982.4585, 37.2243, 5.8297, 1.8295) are given in Table 6. 
The three regression models predicted the aircraft weight with 2 percent error. It was reduced to less than 1 percent 
for the NN technique. For the compressor discharge temperature, the error in regression and NN methods averaged 
0.1 and 2 percent, respectively. The average error in the field length constraints ranged between 1 and 3 percent for 
both approximators. However the error was positive for the regression method, while it was negative for the NN. 

The error in approach velocity was similar to field length constraints. The error was higher for the positivity 
constraints g 4 , g 5 , and g 7 . The solution fidelity was about the same for the small, standard, and the large models. 

To further assess the overall performance of the approximators the errors in the aircraft weight is calculated at 
101 design points for engine thrust (in the range 28 to 35 kip), wing area (1800 to 2200 fe), and turbine inlet 
temperature (2900 to 3100 °R). The mean errors and the standard deviations for the three models is given in Table 7. 


Table 6. Performance of the approximators during analysis 


Response variables 

FLOPS 

solution 

Regression method 
Percent error 

Neural network technique 
Percent error 

Small 

Standard 

Large 

Small 

Medium 

Large 

Aircraft weight in Ibf 

204 725.75 

1.85 

1.67 

1.60 

-0.21 

-0.47 

-0.99 

Approach velocity in 
kn 

112.73 

0.91 

0.83 

0.80 

-1.93 

-0.09 

-2.18 

Takeoff field length in 
ft 

5490.55 

3.66 

3.09 

2.91 

-1.12 

-1.59 

-2.03 

Landing field length 
in ft 

5173.08 

0.94 

0.88 

0.85 

-1.94 

-2.11 

-2.21 

Missed approach 
thrust in Ibf 

3766.80 

-14.48 

-12. 82 

-12.22 

-5.76 

-3.66 

-1.91 

Second segment climb 
thrust in Ibf 

8516.09 

-5.29 

-4.67 

-4.45 

-2.68 

-1.91 

-1.25 

Compressor discharge 
temperature in °R 

1383.64 

-0.07 

0.20 

0.01 

-2.03 

-1.84 

-2.06 

Excess fuel in Ibf 

4237.26 

-77.59 

-70.03 

-66.94 

4.77 

17.08 

24.87 
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Table 7. Percent absolute error in weight over certain design variable ranges 


Variables, range, model 

Regression method 

Neural network technique 

Mean 

Standard deviation 

Mean 

Standard deviation 

Thrust (28 to 35 kip) 





Small 

0.90 

0.15 

1.05 

0.69 

Standard 

1.22 

0.28 

0.92 

0.61 

Large 

1.50 

0.16 

1.01 

0.66 

Turbine inlet temperature (2900 to 3100 °R) 





Small 

0.73 

0.26 

2.08 

1.19 

Standard 

0.81 

0.31 

2.09 

1.21 

Large 

1.10 

0.34 

2.11 

1.21 

Wing area (1800 to 2200 ft2) 





Small 

1.16 

0.22 

1.30 

1.05 

Standard 

1.10 

0.11 

1.12 

0.96 

Large 

1.44 

0.07 

1.16 

0.91 


Both approximators produced about a 1 -percent mean error for the variables, except for a 2-percent error for 
the turbine inlet temperature by the NN technique. The standard deviation in error with the regression method was 
less than 0.3 percent. This was increased to about 1 percent with the NN technique. The error was comparable for 
the small, standard, and large models. 

In the aircraft design optimization, the FLOPS analyzer was replaced by the approximate models without any 
other change. This combined code was run to obtain optimum solution for the aircraft. The combined solution is 
given in Table 8. CPU time to solution is given in Table 5. A convergence graph that shows the aircraft weight 
versus iteration is depicted in figure 2 for the large regression model. From a comparison of this graph with figure 1, 
which used the FLOPS code, we observed the following: 

1 . Design with the approximator required about double the number of iterations than it did with the FLOPS 
code. However, the time to solution was in favor of the approximator: 1.6 CPU seconds for the regression 

method, versus 2031 seconds for the FLOPS code. The NN used 222 seconds. 

2. The convergence pattern contained oscillations for both the regression method and the FLOPS code. The 
amplitudes of the oscillations in the first cascade algorithm were considerably smaller for the regression method, see figure 1 
and figure 2. However, a cascade algorithm was required for the FLOPS code as well as for the regression method. 

3. The approximator exhibited 1 percent error in the optimum weight of the aircraft. For field length and 
approach velocity constraints the error was less than 2 percent. Error was greater for the positivity constraints, which is 
discussed subsequently. 
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{a} Original positivity constraints. 
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(b) Modified positivity constraints. 

Figure 2. Performance of the large model with Regression method with original and modified 
positivity constraints. 
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2.7 Design Sensitivity Analysis 

Design sensitivity was examined for the aircraft to land and takeoff on shorter and longer runways. The length ranged 
between 4500 to 7500 ft with 6000 ft as the nominal value. Other parameters were retained at their nominal values. 

The optimum solutions were depicted in figure 3. Optimum aircraft weight versus the field length obtained by the three 
methods (FLOPS, NN, and regression) was shown in figure 3(a). Likewise the overall pressure ratio and second segment 
climb thrust were given in figures 3(b) and (c), respectively. The approximators exhibited less than 1 percent error in the 
aircraft weight. Aircraft weight was increased for shorter field length and it was decreased for longer length, as expected. 

The NN and regression methods exhibited 0.34 and 0.61 percent error, respectively. Overall pressure ratio (OPR) constraint 
was graphed in a magnified scale in figure 3(b). The approximators hardly exhibit any deviation from the FLOPS solution. 
Observe however, the discontinuity in the OPR constraint. The regression method had a tendency to hug the data points while 
NN exhibited a propensity to follow a mean path. The discontinuity will adversely affect the aircraft optimization when the 
FLOPS code was used. The NN should experience no limitation in design optimization. The performance of the regression 
method for design optimization was expected to be intermediate between the FLOPS code and NN method. Behavior of the 
second segment climb thrust constraint was similar to that for OPR. The regression method closely followed the constraint 
while NN traced an average path. 
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(c) Second segment climb thrust. 


Figure 3. Subsonic aircraft sensitivity analysis with respect 
to field lengths for the small model. 
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2.8 Positivity Constraints 


Three constraints of the problem: missed approach gradient thrust, second segment climb thrust, and excess fuel capacity 
restricted the associated parameter to be positive. The parameters were allowed to approach zero. In a modified case the 
parameters were pushed away from zero, through specified lower bounds of 500 lbf. The optimum solutions for four different 
situations were given in Table 8. From a comparison of the FLOPS modified case to the base design we observed: 

1 . Aircraft weight: The modified regression solution matched the base solution with a 0.58 percent error. The errors 
with the original regression method and FLOPS code solutions were 1 .2 and 0.2 percent, respectively. 

2. Design variables: The modified regression solution for engine thrust and turbine entry temperature exhibited 0.5 and 
0.2 percent error, respectively. The quarter chord sweep angle variable exhibited the most deviation: 16.0 percent for 
the modified regression versus 21 percent for the original solution 

3. Constraints: There was little error in constraints between the modified regression solution and the base design. 
Constraints g4 and gi became active for both the FLOPS and regression methods. For the g5 constraint, the error was 
0.3 percent for the modified regression case. Between the original FLOPS and original regression cases, it was 13 
and 20 percent, respectively. 

3. Overall, the modified regression solution exhibited a closer match with the base design. 


Table 8. Optimum solution with original and modified positivity constraints 



FLOPS 

original 

Regression 

original 

FLOPS 

modified 

Regression 

modified 

Aircraft weight in lbf 

199046.9 

196965.1 

199395.3 

198240.8 

Design variables 





Wing aspect ratio (DV1) 

8.6 

8.8 

8.8 

8.8 

Engine thrust (DV2) 

31408.5 

30179.5 

32346.6 

32181.1 

Wing area (DV3) 

1899.8 

1915.0 

1851.3 

1834.7 

Quarter chord sweep angle (DV4) 

16.0 

16.7 

20.2 

17.0 

Thickness to chord ratio (DV5) 

0.1 

0.1 

0.1 

0.1 

Turbine inlet temperature (DV6) 

3100.0 

3100.0 

3100.0 

3094.5 

Overall pressure ratio (DV7) 

40.5 

38.0 

40.5 

38.0 

Bypass ratio (DV8) 

6.1 

5.0 

6.1 

5.1 

Fan pressure ratio (DV9) 

1.8 

1.8 

1.8 

1.8 

Constraints 





Landing approach velocity (gi) in kn 

119.0 

117.9 

120.7 

120.9 

Takeoff field length (g 2 ) in ft 

6000.0 

6000.0 

5998.3 

6000.0 

Landing field length (g3) in ft 

5479.4 

5425.8 

5562.8 

5573.7 

Missed approach gradient thrust (g4) in lbf 

3746.1 

3136.1 

5000.0 

5000.0 

Second segment climb thrust (gs) in lbf 

8385 

7721 

9619 

9594 

Compressor discharge temperature (g6) in°R 

1429.0 

1405.2 

1428.3 

1405.9 

Excess fuel capacity (g7) in lbf 

0 

2553.8 

500 

500 
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3.0 Eccentric Design 


Design optimization process may have a propensity to converge to an eccentric solution that may be difficult to manufacture. 
This limitation is seldom elaborated in the optimization literature. This section illustrates the occurrence of such solution and 
technique available in CometBoards to overcome the limitation. The subject matter is treated in three subsections, as 
follows. 

1 . Occurrence of eccentric optimum solution 

2. Technique to alleviate eccentric solution 

3. Possible cause of the limitation 

Occurrence of eccentric optimum solution 

Occurrence of eccentric optimum solution is illustrated considering three examples: 

Example 1 — Structural design of an engine component. 

Example 2 — Profile optimization of a cantilever beam. 

Example 3 — Design of a cylindrical shell. 

An example of an indeterminate truss is not included because it is well known that its optimum solution can converge to a 
determinate truss because of the infeasibility of full stress design of indeterminate trusses (ref. 13 and 14). The infeasibility 
theorem is applicable to stress limited design as well as to a design problem with both stress and displacement constraints 
(ref. 15). Consider for example a n-bar truss with r redundant members. For its design under a single load condition areas of 

r-number of members of the truss will reduce to zero (H. = 0, i = 1, 2, r ) . The bars with non-zero areas form the load 

path. An increase in the number of load condition may alleviate the limitation. The infeasibility theorem has not yet been 
extended to flexural structures like continuous beams and frameworks. 

3.1 Example 1: Structural Design of an Engine Component 

An exhaust nozzle of a mixed-flow turbofan engine exhaust nozzle referred to as a “downstream mixing nozzle” shown in 
figure 4(a) is the first example. It is the nozzle for a High Speed Civil Transport aircraft that has to operate at a cruise speed 
of 2.4 Mach in a 5000 nautical miles range. It is fabricated out of rear and forward divergent flaps, rear and forward 
sidewalls, bulkheads, duct extensions, six disk supports, and other components as sketched in the figure. The design 
complexity of the nozzle is increased with flight Mach number, pressure ratio, temperature gradient, dynamic response, and 
degradation of material properties at elevated temperature. The flap is made of Rene 125 material with a Young’s modulus of 
30.4 million psi, a Poisson ratio of 0.3, a density of 0.308 lb/in. 3 and an allowable strength of 1 17 ksi. The flap shown in 
figure 4(b) has a length of 96 in. and a width of 72 in. It is supported by two variable-depth edge beams with a maximum 
depth of 14 in. A grid with a spacing of 12 in. by 12 in. and a depth of 2 in. stiffens the flap. MHOST (ref. 16) and 
MSC/NASTRAN analyzers were used for the static as well as the dynamic analyses of the flap. A finite element model with 
5593 degrees of freedom and 946 QUAD4 elements was used for analysis because both methods produced acceptable values 
for stress, deformation, and frequency responses. For this problem, the thickness of the edge beams, stiffeners, and skin panel 
were considered as the design variables. The flap was designed for minimum weight for von Mises stress, local buckling, 
displacement, and frequency constraints. The problem had 946 stress and 946 buckling constraints, as well as four 
displacements and one- frequency limitations. The design of the flap, which is depicted in Figure 4(c), was obtained using 
three optimization algorithms. All three methods produced the same optimum weight of 1448.2 lbf. At the optimum, the 
frequency at 40 Hz and displacement at 1 in. were the active constraints. Stress and local instability were passive constraints. 

Its animation at the 40 Hz frequency was examined, and one frame is depicted in figure 4(c). The animation exhibited a local 
frequency condition. The edge beams vibrated with significant amplitude, while the response of the rest of the structure was 
insignificant. In other words, only a small part of the structure carried a major portion of the load. In this design, the edge- 
beams became more failure prone than other parts of the structure. The optimum design camouflaged a local solution. 
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[a) Downstream mixing nozzle. 




(c) Fundamental mode animation for 
intermediate solution. [Weight = 1448 lb) 


(d] Fund amenta] mode animation for 
optimum solution. [Weight = 1204.7 lb] 



[b] Rear divergent flap. 


(e) von Mises stress at optimum. 


Figure 4. Design of a rear divergent flap of a downstream mixing nozzle 
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3.2 Example 2: A Cantilever Beam 


The cantilever beam is depicted in figure 5 with two depths. The problem was to generate the optimum profile or depth along 
its length. A total of nine depth variables along its length were considered to calculate its optimum profile. The beam is made 

of aluminum with a Young’s modulus ^i? = 10.1xl0 6 psij, a Poisson’s ratio (v = 0.3) , and a weight density 

(^p — 0 . 1 Ibf / in 3 ^ . It is 240 in. long and 6 in. was the initial depth. The beam weight was the merit function and stress 

and displacement were its behavior constraints. For the purpose of analysis, the beam was modeled with six 20-node- 
hexahedral elements of the IFM/ Analyzers code (ref. 17). The finite element model with 160 displacement degrees-of- 
freedom was used for the analysis model because it adequately predicted the stress and displacement responses of the beam. 
The profile optimization problem had nine depth design variables, and its nine constraints included eight stress and one 
displacement limitation. The profile optimization of the cantilever beam was attempted with uniform upper and lower bounds 
for the design variables at 0.5 in. and 15 in., respectively, see left side of Table 9 . The optimization problem was solved 
using a quadratic programming algorithm. The solution obtained was depicted in figure 6, with solution in Table 9. The 
optimum weight was 1697.5 lbf, and the root stress was the only active constraint. An odd-shaped profile shown in figure 6 
was obtained. The profile was peculiar because the free end, corresponding to a zero stress condition, had a depth of 5.14 in. 
instead of an anticipated lower bound depth of 0.5 in. The optimum solution most likely challenged analysis assumptions, 
and the optimization iterations encountered difficulties. The situation did not improve when a different optimization 
algorithm or a cascade strategy was employed. 
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Figure 5. A simple cantilever beam 
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Optimum depth in inches 



Profile optimization of a cantilever beam for stress and displacement constraints 

Figure 6 


Table 9. Intermediate and final solutions for the beam profile 


Intermediate Solution 

(only stress constraint at the root was active) 

Final Solution 

(Tip displacement and root stress were active 
constraints) 

Design 

variables 

Lower 

bound 

Initial 

design 

Upper 

bound 

Optimum 

design 

Lower 

bound 

Initial 

design 

Upper 

bound 

Optimum design 

1 

0.5 in. 

10.0 

15.0 

14.68 

10.0 

12.0 

13.0 

11.72 

2 

0.5 

10.0 

15.0 

2.98 

9.0 

8.0 

12.0 

9.00 

3 

0.5 

10.0 

15.0 

5.04 

6.0 

7.0 

10.0 

6.22 

4 

0.5 

10.0 

15.0 

6.27 

5.0 

6.0 

9.0 

6.17 

5 

0.5 

10.0 

15.0 

5.41 

4.0 

5.0 

8.0 

5.02 

6 

0.5 

10.0 

15.0 

0.5 

3.0 

4.0 

7.0 

3.67 

7 

0.5 

10.0 

15.0 

2.49 

2.0 

3.0 

6.0 

2.0 

8 

0.5 

10.0 

15.0 

2.67 

1.0 

2.0 

5.0 

1.32 

9 

0.5 

10.0 

15.0 

5.14 

0.5 

1.0 

4.0 

0.50 
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(a) Cylindrical shell. 


Strategy Optimum weight, lb 



Figure 7. A cylindrical shell 


3.3 Example 3: A Cylindrical Shell 

A cylindrical shell with rigid diagrams under two line loads is shown in figure 7 (a). It is made of 
steel with a Young’s modulus = 30 x 10 6 psi) , a Poisson’s ratio (V = 0.3) , and a weight density 

[p = 0.289/6// in 3 ) . It has a radius of 100 in. and length of 200 in. Because of symmetry, only one-eighth of the 

structure was considered for design. A finite element model with 100 QUAD4 elements of the MHOST analyzer was 
considered adequate to predict its response. Its design was cast as an optimization problem with weight as the objective 
function and constraints were imposed on stress and displacement. The thickness of the cylinder along its length and 
circumference were considered as design variables through a profiled depth link factor to provide a heavier design under the 
load. The one eighth shell model was divided into four substructures along its length. The substructures were clustered 
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to obtain three and four subproblems for sequential and parallel computations, respectively. An additional fourth subproblem 
was required to avoid convergence difficulty in parallel calculations. An irregular optimum solution was obtained for the 
shell as shown in figure 7 (b). The depth increased at the center section of the structure while it was reduced towards the 
boundary. The irregular or stepped design represented an eccentric optimum solution. 

3.4 Irregular local optimum solutions 

The three examples, namely the flap (shown in figure 4), the cantilever beam, (shown in figure 6) and the cylindrical shell, 
(shown in figure 7 b) demonstrated that design optimization process can converge to irregular local design solutions. It was 
further observed that the difference may be very little between the optimum weights of the odd design and the regular 
solution. For example, the optimum weight was 1697.5 lbf for both solutions of the cantilever beam, see figure 6. For the 
cylindrical shell, the two optimum weights, for the odd shaped design and the regular design were 1 161.95 lbf and 1 154.1 
lbf respectively. The difference of 0.676 percent in the solutions could be considered negligible especially for the problem 
complexity. The depth differed substantially between the two solutions for the cylindrical shell example. At the crown, the 
optimum depths of 1.322 in. and 2.471 in. obtained for the two designs differed by 53 percent. In summary the optimization 
process can lead to local solutions with little difference in the value of the merit function or weight but the design variables 
can vary substantially. We discuss two important issues: 

• How to bypass local solutions 

• Cause of local solution 

The procedure available in CometBoards and used to bypass local solutions for several optimization problems is discussed 
first. Then we will initiate a theoretical discussion on cause of multiple local solutions. 

3.5 Technique to bypass local optimum solution 

Consider the flap structure shown in figure 4. The local optimum solution of the flap was modified through an examination of 
the animation of a series of designs. The final modified configuration was obtained by increasing the depth of a single edge 
stiffener to 6 in. from its original 3 -in. depth. This configuration was optimized. The increase in the material for this stiffener 
was compensated for by a reduction in the thickness of the edge beam by 58 percent between the two configurations. The 
dynamic animation of the flap and the von Mises stress distribution are shown in figures 4(d) and 4(e), respectively. The 
structure vibrated in a breathing type of mode, or the entire flap responded as a single unit. Local modes were eliminated. 

The optimum design also displayed a full stress condition for a major portion of the flap as shown in figure 4(e). The 
optimum weight at 1204.7 lb. was 20 percent lighter than the original configuration. Animation assisted optimization reduced 
the weight and improved the overall design of the flap. Animation also improved the design of a Spacer Truss of the 
International Space Station. It is not discussed here but elaborated in reference 17. The truss had an open face, which was a 
functional requirement for the on-orbit integration of the spacer truss to other modules of the International Space Station. 
Animation of the optimum design exhibited unacceptable excessive distortion for the open face. The face was braced 
temporarily by a pair of turn-buckles, which weighed 2.5 lbf. during its launch in the space shuttle. The bracing could be 
removed on orbit to allow its integration with the space station. The modified configuration, with tum-buckle bracing, was 
optimized again. It was observed that the open face did not suffer excessive deformations. The final optimum design had a 
weight of 3 16.63 lbf, which was 36 percent lighter than the manual design of 500 lbf. The fundamental mode of vibration of 
the final design was a breathing type of mode in which the entire structure participated as a single unit. In other words, the 
load path was well distributed amongst all members of the spacer structure. Without the use of animation, the generation of a 
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lighter manufacturable optimum structure would have been difficult. The animation and optimization combined tool 
successfully generated an optimum manufacturable design for the spacer truss. Augmentation of animation into the design 
tool CometBoards has successfully bypassed local eccentric solution to generate manufacturable optimum designs for a 
number of industrial structural components. 

The profile optimization of the cantilever beam, shown in figure 7, was attempted with uniform upper and lower bounds for 
the design variables at 0.5 in. and 15 in., respectively. The optimum profile was peculiar because the free end, corresponding 
to a zero stress condition, had a depth of 5.14 in. instead of an anticipated lower bound depth of 0.5 in. The problem was 
solved successfully with manual intervention. Instead of uniform upper or lower bounds and a single optimization algorithm, 
the design was cast as a sequence of subproblems. Upper and lower bounds and initial design through engineering 
intuition were progressively changed for the subproblems. This procedure produced a converged solution that is shown in 
figure 6. It had the same minimum weight of 1697.5 lbf, which was identical to the odd-shaped design. Its root stress and tip 
displacement were the active constraints. For this problem, optimization algorithms converged to two distinct local solutions 
with equal minimum weight. Industry will be more inclined to accept the monotonically profiled beam than the odd- shaped 
design. 

The optimum design of the cylindrical shell problem revealed two local solutions with about the same weight but different 
values for the design variables. At the optimum, the regular optimization and 

the subproblem strategy produced a different number of active constraints. For this problem the attractiveness of subproblem 
strategy was challenged because the design thus obtained was not suitable to manufacture. We are not aware of a scheme to 
alleviate this limitation of the subproblem solution strategy. 

An animation-assisted optimization successfully solved the flap design problem. The beam profile problem was solved 
through an incorporation of engineering intuition. Generating a regular manufacturable design using a subproblem 
optimization scheme still remains a challenge. Bringing optimization methods to their rightful industrial environment from 
the academic plane requires the combined effort of designers and optimization researchers. 

3.6 Cause of Local Optimum Solutions 

It is rather difficult to ascertain ah possible reasons that can lead to local solutions. There exits little literature on the subject. 
We will attempt to illustrate the cause through two simple examples. First example is a three bar truss shown in figure 8 and 
the second example consists of two beams depicted in figure 9. The truss example will discuss strength constraint and 
eccentric design. Stiffness constraint and eccentric design will be discussed considering the beam examples. 
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Figure 8. — Design of indeterminate truss, (a) Determinate, 
(b) Indeterminate. Dimensions are in inches. 



(a) Determinate cantilever beam (b) Indeterminate clamped beam 

Figure 9: Sensitivity and stiffness constrainst 


1) The Three Bar Truss 

Traditionally the concept developed to design a determinate structure is also used to design an indeterminate structure. 
Consider for example a determinate truss. The concept of fully stressed design (FSD) works very well for a determinate truss 
shown in figure 8 (a). Extending the FSD concept to design an indeterminate truss can become problematic (see ref. 13-15) 
because an indeterminate truss cannot be fully stressed. Rules to design a determinate structure should be modified prior to 
their application to design an indeterminate structure. The compatibility conditions are key ingredients that should be 
accounted for the modification because the CCs imposed implicit constraints in the design of an indeterminate structure. The 
implicit constraints are intrinsic to the design of an indeterminate structure and are independent of the choice of an analysis 
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method used in the design calculations. The nature of the constraints is easily studied via the IFM, however. The rules thus 
generated should be adopted even when the stiffness method analysis tool is used in design calculations. Cilley (ref. 18), in 
1900 discussed some aspects of design of an indeterminate truss. Reinschmidt et al. (ref. 19) attempted an analytical 
explanation through an illustrative example. Gallagher and Zienkiewicz (ref. 20) has also touched upon the issue. Some of 
our observations on the topic are reported in references 14, 15, and 21. It turns out the CC has considerable influence on the 
fully stressed and fully utilized designs as well as in the optimization methods. The central design issue is illustrated by 
considering the truss example. Extension to indeterminate beams, framework, and other types of structures is straightforward, 
though it would require extensive algebraic manipulation. Both traditional design and optimization concepts are examined for 
an indeterminate truss. The discussion is separated into two topics: 

• Design for strength 

• Design for stiffness 

Design for Strength Design for strength attempts a full stress state for all members of the truss. Consider a truss with n bars. 
For all bars in the truss a full stress condition is defined as 

(T = CT 0 (i = 1,2, ...,n) (6) 


where o is the stress in the i th bar and <T 0 is the strength. A stress ratio rule yields the area Ai of the i th truss bar as a ratio of 


r 


member force F, to strength: 
indeterminate truss. 


4 


. The rule works well for a determinate truss but becomes problematic for an 


Design of a determinate truss . — Consider a determinate truss with n bars and the same number of bar areas. The bar areas can 

be calculated from the governing equation of the IFM^^jjT 7 } = | j. by expressing it in areas A t = FJa 0 to obtain 

[iS ] {^4 } = {P }. For a determinate structure the governing matrix is equal to the equilibrium matrix ([, S] = [B]), which is 
independent of bar areas. In other words, a single inversion without iteration yields the areas: 


Ml = abs{[5'*]" 1 {P*}} 


( 7 ) 


The formula in equation (7) can be extended to include buckling limitations. 

Consider the two-bar determinate truss shown in figure 8a. It is made of steel and is subjected to two load components, P x = 
10 kip and P y = -20 kip. The bar force is linked to the cross-sectional area A through the yield stress as F\ = <JqA\ and F 2 = 
OqA 2 . The IFM governing equation degenerates into the equilibrium equations. The EEs are modified and solved to obtain the 
design. 
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MM = M =• MM=( p l =• 

or {^} = abs|W[5]" 1 {P}| (9) 

ll CT oJ J 

The bar areas are determined from equation (9) in a single step because the equilibrium matrix [B] is independent 
of the design variables. The numerical values for the areas are as follows: 



weight = 43.4 lbf (10) 


Design of an indeterminate truss . — The truss is made indeterminate by adding a third bar as shown in figure 
8(b). The following design equation is obtained when the determinate design concept is extended to the 
indeterminate truss: 


[S]{F} = {/>*} 

or, {4 = abs|^j[5] _1 {p*}| 


( 11 ) 


Equation (11) has to be solved iteratively since the [S] matrix is a function of bar areas because of the flexibility 
matrix [G] in the compatibility condition [C][G]{F} = {0}. Iterative solution of equation (11) yielded the 
following design: 

0.7f 

0.50 >in. 2 and weight = 43.4 lbf 

0.00 

An extension of the determinate design concept to an indeterminate structure returned with a zero area for the 
third bar. Furthermore, design and weight are identical for both determinate and indeterminate structures. This 
observation is not coincidental but is common to all indeterminate trusses. The implicit compatibility constraint 
differentiated between the designs of the determinate and indeterminate structures. This implicit compatibility 
constraint in deformation, force, and stress variables {(3), {F}, and {a}, respectively, is as follows: 
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Pi - V2P2 + P3 - 0 


(12 a) 



Gl-a 2 +a 3 =0 (12c) 


The compatibility constraint is violated (ref. 14) when a fully stressed state (gi = a 2 = a 3 = a 0 ) is attempted: Gi - 
a 2 + a 3 =^> 20 - 20 + 20 ^ 0. Because of the violation, a bar area of zero is obtained for the third bar (A 3 = 0) of the 
determinate truss. This observation can be extended for a general indeterminate truss with n bars and r redundant 
members. For such a truss ( n , m) there are m = n - r equilibrium equations and r = n-m compatibility conditions. 
A fully stressed condition cannot be achieved for an indeterminate truss ( n , m). An attempt to fully stress this truss 
would degenerate it to a determinate truss with areas of zero for r number of bars: Aj = 0 where j = = 

3.7 Eccentric optimum design and the compatibility condition 

The indeterminate truss with r-number of bar areas, either equal to zero or the minimum bounds for areas is an 
eccentric optimum structure. Compatibility condition, which became the implicit constraints, become the cause 
of eccentric design. The impact of compatibility condition has not yet been studied in design of flexural members 
like beam and frameworks. It is believed that hinge formation may take place in the optimum design of 
continuous indeterminate beams and such determinate structure may also be considered as the eccentric optimum 
solution. 

Linking of design variables . — The situation associated with bar areas of zero in the design of an indeterminate 
truss can be alleviated by linking the design variables. The number of independent variables must be reduced to m 
for an n-bar truss with r redundant members. For the three-bar truss (see figure 8(b)), we discuss two of several 
choices to link areas: link the areas of bar 1 and bar 3 (A u A 2 , A 3 = A\); likewise link A 2 and A 3 . The designs 
obtained for the two cases are as follows: 

0.64 

{Aj A 1 and ^ 3 =< 0.66 > in. 2 and weight = 71.5 lbf 
0- 64 

0.50 

{A} A 2 and ^ 3 = < 0.71 > in. 2 and weight = 63.8 lbf 
0.7! 
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The weight of the structure is increased to 71.5 or 63.8 lbf from the weight of 43.4 lbf for the determinate system, 
when areas are linked. Furthermore, the weight of the structure depends on the choice of the link strategy. For the 
three-bar truss, weight was increased by 65 percent from the determinate truss when areas A\ and A 3 were linked. 
The increase was 47 percent for the link condition A 2 = A 3 . 

The link strategy has considerable influence on the weight of the structure. However, it is not straightforward to 
determine the best link strategy. The bandwidth of the CC provides some analytical guidance to develop a link 
strategy. Consider for example the 20-bay truss shown in figure 10. This truss has 101 members with 21 
redundancies and 21 CCs. The bandwidth of the CCs are 6 for 20 of them and 20 for 1 of them. The bandwidth of 
six pertains to a bay in the truss; that is, the first bay has six bars (1 to 6). The areas of these six bars are 
candidates for linking. A simple strategy is to link the diagonal bars (A 3 = A 4 ). Likewise link bars (A$ = A 9 , A 13 = 
A i 4 ,. . ., A 9 % = A 99 ) for other bays. The remaining CC pertains to the 20 bars of the bottom chord. A simple strategy 
is to link all such bars ( A 2 = A 7 = A n = ... A 97 ). The compatibility bandwidth provided a guideline, but there is no 
unique strategy. It may be tempting to employ linear programming for the purpose. The stiffness method cannot 
provide a link strategy because of the nonavailability of the bandwidth of the CCs. The CC explained why an 
indeterminate structure cannot be fully stressed. The bandwidth of the CC can be used to develop a linking 
strategy. 



Design for Stiffness: Design for stiffness included displacement limitations in addition to the strength constraints. The 
design obtained for strength remained adequate if displacement became a passive constraint. The design has to be modified 
when displacement constraints become infeasible. The design can be prorated uniformly to satisfy the displacement 
limitations. The prorated technique can produce an overdesign condition because of an underutilization of the strength. The 
following discussion is an attempt to alleviate the overdesign condition. 

Stress-displacement relationship . — Displacement and stress are dependent upon one another, regardless of the bar area of a 
tmss. The dependency relationship (ref. 15) is derived for the three-bar tmss shown in figure 8b. 

{PH^} =[*]>} (13a) 
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(13b) 


(13c) 

"3 -£(*+*) 

(13d) 

II 

1 

fra 1 a 

i 

Q 

(13c) 

*2 =-§(*2) 

(I3f) 


The displacement and stress relationship can be generalized in the following formula: 

Pa 

°k = X C jk X j ( 14 ) 

j=P\ 

The coefficients c jk are independent of the design variables, here bar areas. The actual number of displacement 
components Xj in equation (14) depends on the column bandwidth of the equilibrium matrix [ B ] that is a small 
number. The maximum bandwidth is four for a two-dimensional truss. It is six for a space truss. In other words, a 
small number of displacements are related to a stress component. For the three-bar truss, displacement depended 
on one or two bar stresses. 

One approach is to modify the stress allowable or strength for selected members using equation (14). Then 
proceed with the fully stressed design. Let us assume the design for strength obtained for the three-bar truss 
indicated a 15 percent violation on the limitation specified on the displacement component X 2 . From equation 
(13f), we can satisfy the displacement constraint if the corresponding allowable stress is reduced by 15 percent 

(a 2 o = 0.85 g 0 ). This in turn would increase the area of second bar by 17.6 percent (a 2 qw = 1.176A 2 ld ) • The 

design update requires no iteration for a determinate truss. Iteration may be necessary for an indeterminate truss. 
The rule is easily generalized for a truss with k violated displacements. 

Design formula for displacement limitation . — A general formula to design for displacement limitation can be 
obtained by modifying the displacement-force relationship of the IFM. 

{X} = [J][G]{F} (15) 
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The rectangular matrix \J], is replaced by \_[S] 1 J with appropriate modification to obtain a the displacement 
formula as follows: 

|^| = [[^r 1 ] r [G]{F} = [5 r ]{4 (16) 

Equation (16) is manipulated to obtain an expression for the bar areas for displacement limitations as 
[s]U} = j^2j, or UMsThx} (17) 

Equation (17) is combined with equation (1 1) to obtain formulae for stress and displacement limitations: 


{^} stress =absM 1 V>} 

(18a) 



U} dtap =abs([sr 1 {JP}) 

(18b) 

{4 = max({4 stress ,{4 disp ) 

(18c) 


The formulae in equation (18) are modified for computational efficiency and superior convergence characteristics 
to obtain a method called modified fully utilized design (MFUD). Discussion of the MFUD technique is given in 
reference 13 and is not repeated here. However, results for a set of nine problems are presented in table 10. Each 
problem is solved twice: first by a design optimization method and then by MFUD. 

Consider for example the design optimization of the 60-bar trussed ring. It was subjected to three load cases. It 
had 60 stress and 6 displacement constraints for each load case. The design problem was solved using MFUD and 
an optimization method using sequential unconstrained minimization technique (SUMT). The 60 bar areas were 
linked to obtain 25 design variables. For the MFUD method and SUMT algorithm, the optimum weights were 
308.1 and 309 lbf, respectively. The active constraint set for MFUD method included 18 stresses and 
1 displacement, against 1 1 stresses and 1 displacement for SUMT. The MFUD method generated a good design 
solution with few calculations. The other eight examples in table 10 follow the pattern of the ring problem with 
minor variations. 
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TABLE 10.— NORMALIZED WEIGHT OBTAINED BY SUMT AND MFUD METHODS 


Problem 

Normalized weight 
(with respect to SUMT solution) 

Number of active constraints 

SUMT 

MFUD 

Stress 

Displacement 

SUMT 

MFUD 

SUMT 

MFUD 

3 -bar truss 

1.0 

1.00 

2 

2 

1 

1 

5 -bar truss 

1.0 

1.00 

- 

- 

1 

1 

Tapered 5 -bar truss 

1.0 

1.00 

2 

3 

1 

2 

10-bar truss 

1.0 

1.02 

1 

- 

2 

2 

Tapered 10-bar truss 

1.0 

1.00 

5 

3 

2 

2 

25-bar truss 

1.0 

1.00 

2 

4 

4 

4 

20-bay truss 

1.0 

1.02 

13 

13 

2 

4 

60-bar truss ring 

1.0 

1.00 

18 

18 

1 

1 

Geodesic dome 

1.0 

1.01 

46 

46 

1 

1 


2(a) The Cantilever Beam 

Design for stiffness consideration is illustrated considering the two beams shown in figure 9. The cantilever beam shown 
figure 9(a) was made of aluminum material. The determinate structure was subjected to a tip load of 20 kip. The beam depths 
(di and d 2 ) were the two design variables while the passive thickness variable (f = 4.0 in.) was uniform across its length. 
Other parameters were shown in the figure. The design objective was to minimize weight under strength and stiffness 
constraints. The indeterminate clamped beam shown in figure 9(b) was also made of aluminum with properties identical to 
the cantilever beam. It had a length of 4a = 1 92 in. The uniform thickness of the beam was set to (f = 2 in.). The beam was 
subjected to three transverse load components as shown in the figure with a total magnitude of (P = 1 00 /c/p) . Allowable 
stress limit was (cr 0 =20/cs/)and tip displacement was not to exceed one quarter of one inch(£ <0.25). The depths 
(d v d 2 ) were the two design variables. 

Fully stressed concept yield a design for strength limitation. For the determinate beam in figure 9, the two depths were 
calculated from the induced stresses at locations 1 and 2 as follows: 



My ' 

I )i 


(19) 


where, g 0 was the stress allowable, induced stress was a, the bending moment was M, moment of inertia was /, and y 
represented the half-depth of the beam. Location 1 is at the clamped boundary, while location 2 was at the transition, see 
figure 9(a). Solution to equation 19 yielded the traditional design for the determinate beam. Its weight was back calculated 
and it could be close to the minimum weight even though it was not explicitly used in the design calculation. For an 
indeterminate structure, the equation set 19 have to be solved iteratively. The implicit compatibility constraints may 
produce an eccentric design. This is an area for further research because the effect of the implicit constraints in the 
design of flexural structure have not yet been examined. 
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Local stress constraint 


It is important to observe that stress is a local constraint. The design variable, here depth ( d x ) and the associated response 
variable ( <J X or M x ) occur at about the same location. Stress (cr, ) can be reduced by increasing the depth (^ ) . Likewise 

Stress (cTj ) can be increased by reducing the depth ( d x ) . Proportioning for stress constraints by the stress ratio method is 

very efficient because of the local nature of the stress constraint. Sensitivity of stress constraint for an indeterminate structure 
has two components: a local factor and a system or global factor. The local factor dominates and the second factor can be 
ignored without any consequence (ref.21). In other words a non-gradient based method will perform well for design 
optimization for stress limited design for the determinate beam. This is the reason for the great performance of the stress ratio 
method of structural design for determinate systems. 


Global nature of displacement constraint 


Displacement is a global constraint. It cannot be modified efficiently by changing associated local design. For example the 
tip displacement of the cantilever beam cannot be changed efficiently by increasing the local depth (d 2 ). In other words, it is 
not straightforward to design for a violated displacement constraint even for a determinate constraint. The difficulty is 
illustrated considering two different approaches: an engineering method and a proration technique. 

For the tapered beam, the displacement (£) at the free end can be expressed in closed form as 


S = 


Pa 3 'j 

7 

_i 

E J 

d 3 

V u i 

d 3 

u 2 y 


( 20 ) 


The displacement peaks at the free end which is associated with the design variable ( d 2 ) . It would be inefficient to reduce 
the displacement (£) by increasing the depth (J 2 ) . Likewise it would not be a good design to accommodate the stiffness 

constraint by increasing the depth ) near the support. Both the design variables , d 2 ) have to be changed to obtain 

an efficient design. The experience based ‘Engineering Method’ attempts to satisfy the violated displacement constraint 
intuitively. In this method, all depths are changed uniformly by updating the depths by a factor (k) as: 

dr = k d s x 

d n 2 ew = k d 2 ( 21 ) 

The design for strength constraints only is designated with the superscript ‘s’ in equation 21 . The design with superscript 
‘new’ satisfies both strength and stiffness constraints, assuming k to be greater than unity (k >1.0) . Solution of equations 
(20 and 21 ) for the required displacement limit (S max ) yields the factor (k) as: 
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k = a 


( 22 ) 


f ( p ^ 

7 1 Y| 

V ^max^ y 

V 

1 

m K)’JJ 


Equation 22 can be modified and solved iteratively for a situation with multiple violated displacement constraints. Pick the 
most violated displacement limitation and make it feasible using the procedure discussed. Reevaluate the design to determine 
its feasibility. Satisfy a next violated displacement limitation if any. A feasible design is generated in a few iterations. For a 
general type of structure, a solution to equation 22 or its variant is quite simple in a computer. Engineering method may 
exhibit a tendency to become heavy because the design update is uniform across all members. A little more satisfactory 
design for the problem may be obtained by updating the root depth (c/ 1 ) without any change to (c/ 2 ) . This concept of 
selective updating of the design is discussed next 

Proration Technique: The proration technique is an attempt to alleviate the limitation of the engineering method, which 
may have a tendency to produce an over-design condition. Instead of uniform proration, a weighted proration can be used. 
The weight factors are calculated from the sensitivity of the displacement constraint. It is important to observe that sensitivity 
is required in the design calculation with stiffness or displacement constraints. For the tapered beam the sensitivity can be 
calculated in closed form as: 


ir {s)= - 

5c/ 1 


dd. 


-(S): 


Pa 3 \ 

21 1 

E J 

U 4 J 

Pa 3 ) 

3 1 

E J 

U 4 J 


{0.36, normalized to pr \ = 1 .75} 
{0.19, normalized to pr 2 = 1 .00} 


(23) 


Then equation 23 can be rewritten as: 

d™ w = (pr, = 1 .75) A d* 

d n 2 ew = (pr 2 =1.00) kc/ 2 s (24) 

The weight factors bias the design towards (of 1 ) which is the depth of the cantilever near the support and it is intuitively 

correct. The proration method is discussed in reference 26 for indeterminate trusses. The calculation of sensitivity has been 
simplified and it does not require intensive computation (ref. 21). The engineering method and the proration technique were 
developed considering the example of a determinate structure. The formulations, however, can be used for indeterminate 
flexural structures. Design iterations are required for strength and stiffness limited design for an indeterminate structure. The 
infeasibility of a fully stressed design (ref. 15) has to be addressed when the formulation is applied to indeterminate trusses. 
This design deficiency may not be pronounced for indeterminate flexural structure 

A summary of the designs generated by the methods was given in Table 1 1. For the strength-limited design, the traditional 
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stress ratio method yielded the two depths, that were calculated from moments at the base and at the transition point for the 
nominal strength of (<r 0 = 20 ksi). The optimization methods also converge to the same design point. The design 

(d, = 14.7/n., d 2 = 10.4/n. and weight = 722.9 Ibf ) with a tip displacement of (S max = 2.308 in.) violated the stiffness 
limitation of one inch ( <7 < 1) . 


Table 11. Design of tapered beam by different methods 


Method 

Design variables in 
inch 

Weight 

(ibf) 

Max stress 

Max disp 

d, 

d 2 



Optimality Criteria 

20.4 

12.5 

946.9 

13.8 

1.0 (active) 

Proration Technique 

19.4 

13.8 

955.4 

11.4 

1.0 (active) 

Engineering Method 

25.0 

10.4 

1018.1 

20.0 

1 .0 (active) 

NLPQ Algorithm 

20.4 

12.6 

951.9 

13.6 

1.0 (active) 

Cascade Strategy 

20.5 

12.5 

951.9 

13.8 

1.0 (active) 


The design had to be adjusted for the feasibility of the displacement limitation. Generation of a design for both stress and 
stiffness constraints was not straight forward even for the simple determinate structure. The maximum displacement (<7 max ) , 
which occurs at the tip of the cantilever, can be expressed in terms of the two depths as: 




7 46.5 


r 7_ 


(25) 


The displacement function was graphed in figure 11. The x-axis represented a pair of (c/ 1 and d 2 ) or the design for which 

(<7 max = 1/n.) . It was observed that there were an infinite number of designs that satisfied the displacement limitation. 

Design generated by different methods are listed in table 1 1 and marked in figure 1 1. For the problem, the optimality criteria 
method produced a design with the least weight of 946.9 lbf. The weight calculated by the prorated design method at 955.4 
lbf was about one percent heavier than the solution by the optimality criteria approach. The design was heavier by 7.5 percent 
by the engineering method. The quadratic programming technique (NLPQ), and the cascade strategy converged to design 
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with a optimum weight of 95 1 .9 lb which is less than one percent heavier than the optimality criteria solution. A lighter 
design was produced in the proration method that changed both the depths. Design became heavier when only the depth of 
the support member was changed. In other words, intuition was not satisfactory for the determinate problem. In summary, the 
maximum variation in the weight at 7.5 percent was considered acceptable. Deterministic design methods produced (29 and 
33) percent variation in the design variables (di and d 2 ), respectively. In other words, the variation in the depths was 
substantial for the simple problem. For all design points, the displacement was the active constraint, while the two stress 
constraints were passive. 
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Design variables (di,d 2 ) in inch 


Figure 11. Design for a tapered beam for strength and stiffness limitations 


2(b) Clamped Beam 


The indeterminate clamped beam shown in figure 9(b) was made of aluminum with properties identical to the tapered beam 
discussed earlier. It had a length of 4a = 1 92 in. The uniform thickness of the beam was set to (f = 2 in.) . The beam was 
subjected to three transverse load components as shown in the figure with a total magnitude of(P = 100 kip) . Allowable 
stress limit was (cr 0 =20 ksi ) and tip displacement was not to exceed one quarter of one inch (S < 0.25) . The depths 
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( d v d 2 ) were the two design variables. 


The stress ratio method produced the design for strength: (c/ 1 = 20. 2 in., d 2 = 11 .5 in.) with a weight of 608.4 lbf. The design 
converged in five iterations because of indeterminacy. The maximum displacement was (8 = 0.38 in.) but the requirement 

was(£ < 0.25) . Solutions generated by different methods for both stress and displacement constraints were listed in table 12. 
The proration technique converged to a weight of 699.1 lbf. The weight by engineering method was 707.5 lbf. The weight 
calculated by the optimality criteria method was intermediate to the proration and the engineering method. The design 
generated by the optimization methods converged to a structure with four percent higher weight than that was calculated by 
the traditional methods. The variation was about fifteen percent for the design variables. The formula for the transverse 
displacement at the center of the beam can be written as: 


8 =138.24 

max 


2d & 2 +21c/ 1 3 d 2 3 + d° 


(26) 


d'tdKd't +dl) 

The displacement constraint (b' max = 0.25 in.) against the beam weight could be graphed to obtain a figure 


that would look similar to figure 11 and it was not repeated. The message however is the same. There are 
an infinite number of designs that satisfy the displacement constraint ( g = J max - 0.25 = 0) 


Table 12. Design for clamped beam by different methods. 


Method 

Design variables in 
inch 

Weight 

(lbf) 

Max stress 
(ksi) 

Max disp in inch 

d i 

d 2 

Proration Technique 

23.5 

12.9 

699.1 

15.13 

0.25 

Engineering Method 

21.5 

15.3 

707.1 

19.0 

0.25 

Optimality Criteria 

23.8 

12.8 

702.7 

15.4 

0.25 

NLPQ 

24.5 

13.5 

728.6 

14.2 

0.25 

Cascade Strategy 

24.5 

13.5 

728.6 

14.2 

0.25 
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4.0 Design sensitivity and singularity in structural optimization 

Design sensitivity is central to most optimization methods. It can be a major contributor to the number of calculations in 
optimization. The computation of efficient design sensitivity for structural problems has drawn considerable attention. NASA 
organized a conference on the subject matter (ref. 23). It is also well documented in the literature (23-256). General-purpose 
codes provide for the calculation of sensitivity for stress and displacement constraints. The automatic differentiation of the 
FORTRAN code, ADIFOR has also been suggested to calculate sensitivity. In such a circumstance, we ask and attempt to 
answer a question about the precision of the sensitivity of stress and displacement constraints in the design optimization of an 
indeterminate structure: “Is the optimization process robust when the design sensitivity matrix is highly accurate?” The 
contrary may be true. The performance of an optimization method can be improved when the analytical sensitivity is replaced 
by a determinate approximation (ref. 26). The approximate sensitivity matrix is not only adequate, but it should be preferred 
in the design calculation of an indeterminate structure. Optimization, in other words, requires sensitivity, but approximate 
gradients are quite satisfactory. To illustrate the benefits that accrue from the approximations, the authors used several 
indeterminate trusses as numerical examples because design optimizations have been completed for such structures. The 
concept, however, should be extendable to other types of structures, such as beams, framework, and shell structures. 

Consider a truss that is made of n bars with r dependent members. The n bar areas are treated as the design variables. 
Approximate sensitivity works well because of three attributes special to an indeterminate truss. 

(1) Dependent stresses: In an indeterminate truss, r out of n bar stresses { } are dependent. Stresses are dependent 
because of the r compatibility conditions, which can be written as |^C^{a} = {0}. The r x n sparse matrix [^Cj is 
independent of the design variables. 

(2) Dependency of stresses and displacements: The m = n-r number of displacements {X} are dependent on the n bar 

stresses: {X} = [ 5 ] {a} . The m x n sparse matrix J is independent of the design variables. 

(3) Active constraints: In structural design, the number of active constraints can exceed the number of design variables. 

In an optimization algorithm, the singularity condition can be alleviated by restricting the number of active constraints to not 
exceed the number of design variables. 

In an optimization algorithm, the calculation of the search direction {d} requires the constraint gradient matrix [Vg]. An 
example of the use of the sensitivity matrix to generate the direction follows: 

M=[e]{v/( 

_ [-f] - [vg][[vg] r [vg]]~‘ [vg] r 

0.5^{[ff]{v/}} 7 '{[ff]{ v /}} 

(27) 


Here/is the objective function, and 


[tf]=W-[v4[?*f [ v dl Vsf 
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The matrix 




becomes singular for each of the three special attributes, items (1) to (3), stated above. An 


optimization algorithm may yield a solution despite the singularity condition because the \Q\ matrix is approximated most 

often; it is seldom calculated in closed form. It is reinitialized into an identity matrix when corruption is suspected. The 
proposed approximate sensitivity of stress and displacement constraints alleviates the singularity condition in the design 
optimization of an indeterminate truss. The solution is reached with fewer calculations because the optimization process 
becomes more robust, and the sensitivity is generated with a trivial amount of computations. The benefit that accrues from 
the use of approximate sensitivity is shown through the solution of a set of problems that were selected from the literature. 
Each problem was solved twice in a controlled environment, first using the closed- form gradient, then with an approximation. 
A comparison of the two optimum solutions quantified the benefit. The underlying cause of the benefit was investigated 
through a discussion of the nature of structural design optimization problems. 


Consider the / h stress constraint. Its closed-form gradient can be expressed as the sum of two factors: 




determinate r \ indeterminate ) 

+ K! 


( 28 ) 


The gradient expression given by equation (28) has to be adjusted for the absolute value in the constraint, which however, 
poses no limitation to the discussion here. The first factor {V ^.} determinate i s applicable to determinate as well as indeterminate 
trusses. This vector has only one nonzero entry, which is the negative ratio of the member force to the square of the bar area 
( -FI A 2 ). The second term {V ^.| mdeterminate accounts for the effect of indeterminacy. It is not a negligible factor. It can be fully 
populated, and its calculation is computationally intensive. The proposition is to drop the second term {V ^| mdeterminate [ n 
design optimization even when it is nontrivial. The nature of the gradient of the displacement constraint is quite similar to 
that of the stress constraint. Again, the proposition is to retain only the simple determinate factor. 

The gradient-approximation concept is illustrated considering torsion of a shaft as an example. The shear stress (r), at a 
distance T from the neutral axis, with an induced torque ( T ) and polar moment of inertia ( j) can be written as: 



Consider the polar moment of inertia J as the design variable. For a determinate shaft the induced torque is independent of 
the polar moment of inertia and gradient with respect to J can be written as: 


y T det ermmate = = _A jy (2%) 

8J J 2 

The gradient matrix for a multiple variables shaft problem becomes a diagonal matrix with elements given by equation (29b). 
For an indeterminate shaft the induced torque T(J) becomes a function of the polar moment of inertia. Equation 29b has to be 
modified as: 
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Vr = 


Deter min ate 


indeterminate 


Determinate 


-I TF 

+ 

r 

VT * 

—jTF 

J 2 


7 


J 2 


(29c) 


The proposition is to retain only the determinate factor, which is easy to calculate and drop the indeterminate factor that is 
computation intensive because for indeterminate structures the torque becomes an implicit function in the design variables. 


An analysis tool is required to calculate the stress and displacement constraints and their sensitivities. Here, the integrated 
force method (IFM) [14] is employed. The structure of the IFM equation is suitable to calculate the closed- form sensitivities 
because the sizing design variables of a structure (here bar areas) are retained in a pristine state in the concatenated flexibility 
matrix [G]. In addition, IFM has two distinct sets of equations. The first set, with internal force as the primary unknown, is 
differentiated to obtain the sensitivity of stress. Likewise, the sensitivity of displacement is recovered by differentiating the 
second set of equations. The adjustment that may be required for the stiffness method of analysis is also discussed in this 
paper. The IFM equations to calculate forces and back-calculate displacements are as follows: 

Internal forces {F 1 } are calculated from the governing IFM equation: 

[S]{F}={P} (30a) 

Displacements {X} are back-calculated from the forces: 

{X} = [J][G]{F} (30b) 


The sensitivity matrix for the stress and displacement constraints for an ft-bar truss with r dependent members is obtained by 
differentiating the IFM equations. The closed- form sensitivity matrix for stress has the following form: 


[Vo] = [Vo 1 ,Vc 2 ,..,Vo )! ] 


T 

- determinate 

r i 

r- t 

- indeterminate 


1 

+ 

E--i 

'• 1 
A--. 


(31) 


The expression given by Equation (31) should be adjusted for the allowable strength prior to its use in design optimization. 
The rank of the n x n sensitivity matrix [V ] in equation (31) is reduced to m = n - r when both terms are retained. The 
recommendation is to use only the first term in equation (31), which is superscripted “determinate.” It is a diagonal matrix 
with full rank n. The proposition is to drop the second term that is superscripted “indeterminate.” The calculation of the 
second term is computation intensive. The closed-form displacement sensitivity follows: 
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[vx] = 





n determinate 


~\T 


+ [[/][G][D]] lnde ‘™ lna,e 


(32) 


The proposition is to use the first factor with superscript “determinate" in design optimization. The 
definitions of the symbols in Equations (31) and (32) are as follows: 


phst 1 



[c] 


The elements of the diagonal matrix 




are given by 


Sii F i 44 


l c \ _ Aii r i _ 

\ /// — T — — 


'• 1 

’• F 


= 


r • -i 




■G-. 

, and 

'F-. 

7 -. : 

A 2 \ 




- ' -1 


4 4 4 

1 are diagonal matrices of dimension nxn whose 


IF- t- 

elements are — , — l — , — — , and F t , respectively. 
A At A; E; 


l l l 


The determinate factor in the displacement sensitivity can be specialized for an ?z-bar indeterminate truss as 


r- 


[ V X jdeterminate = 


-Ft 

A 2 E 


[4 


(33) 




The displacement sensitivity given by Equation (33) is a function of the bar length Young’s modulus E , and areas A 
because displacement is a global variable. Calculation of the determinate sensitivity for the displacement essentially requires 
a back-substitution step with the factored form of the inverse of the S matrix. It is important to observe the similarities and 
differences in the sensitivity expressions of the stress and displacement. 
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1. Sensitivities of both the stress and displacement contain the member forces ji 7 } and the square of bar areas j A} . 

Even for an active displacement constraint, the design in essence is modified through the member force. 

2. The geometry of the truss is not explicitly contained in the sensitivity expression for stress because it is a local 
variable. 

3. Because displacement is a global variable, its sensitivity expression explicitly contains the geometrical or 
configuration parameters, the material property, and the design variables. 


The consequences of using approximate sensitivity in design optimization are demonstrated through the solution of a set of 
problems. Each problem was solved twice, first using the determinate sensitivity, then with the full closed-form expression. 
Problems were solved in a controlled environment on an SGI workstation running the IRIX 6.5 operating system (Silicon 
Graphics, Inc., Mountain View, CA). Identical convergence and stop criteria were used for the optimization algorithm. For 
large problems, we reduced the number of design variables and constraints by utilizing the design variable linking and 
constraint formulation features available in CometBoards. A sequential quadratic programming algorithm, SQP [12], was the 
primary optimizer. This algorithm was supplemented, when required, by a modified method of feasible directions (mFD) and 
a sequential unconstrained minimization technique (SUMT). The examples and solutions given in references ???, are not 
repeated here. 

4.1 Computational efficiency 

The set of examples, (see table 10) converged to the correct solution with the approximate sensitivity. To examine the 
computational efficiency when the approximate sensitivity was used, we solved the three larger problems: the ring, the wing, 
and the 20-bay truss in a controlled environment using the SQP algorithm in an SGI workstation with an IRIX 6.5 operating 
system. The CPU time to optimum solution was measured for both the determinate and the closed- form analytical 
sensitivities. The CPU time to solution is depicted in Table VII. For the 60-bar trussed ring, the optimum solution was 
reached in 2.3 CPU sec with approximate sensitivity. The time increased threefold for analytical sensitivity. For the forward- 
swept wing, the time factor in favor of the approximation was almost eightfold; 8.7 CPU sec with approximation against 67.5 
sec for analytical sensitivity. The time ratio was 6 for the 20-bay truss; 2.5 and 14.8 CPU sec with approximate and analytical 
sensitivities, respectively. Overall, the time to solution increased from threefold to eightfold for the analytical sensitivity. 
Approximate sensitivity increased computational efficiency by several orders of magnitude. 

Table 13. CPU time to solution for large problems with the SQP algorithm 



60 bar ring 

Forward swept wing 

20 bay truss 

Type of sensitivity 

CPU 
time in 
seconds 

2.3 

8.7 

2.5 

Determinate 

CPU 
time in 
seconds 

7.4 

67.5 

14.8 

Analytical 


4.2 Convergence pattern 

The convergence of weight versus CPU time to solution with analytical and approximate sensitivities for the three large 
problems is depicted in figure 11. For each problem, optimization was begun with the same initial design. Both methods 
(with determinate sensitivity and analytical sensitivity) produced similar optimum solutions. Consider the graph in figure 
11(a) for the 20-bay truss. The convergence patterns, with and without approximation, portray undulations that are quite 
similar. However, the convergence is very rapid with the determinate sensitivity. The convergence characteristic is similar for 
the 60-bar trussed ring and the forward-swept wing shown in figures 1 1 (b) and (c), respectively. 


NAS A/CR— 20 13-217748 


151 






n 


c 

O) 

<D 


§ 



0 


(c) 


T 


— Analytical sensitivity 

— Approximate sensitivity 


t 1 1 r 


0 10 20 30 40 50 60 

CPU time, sec 


n 

70 


Figure 11. CPU solution time for SQP 
algorithm, (a) Twenty-bay truss, (b) Sixty- 
bar trussed ring, (c) Forward-swept wing. 
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4.3 Angle between search directions 


The iterative optimization process moved along a search direction ydj from one design point to another. The search 

directions, generated from the gradients, differed for the analytical [d anl jand determinate sensitivities ^J determinate j . If the 

difference between d anl and ^ determinate was small, then the contribution to sensitivity from the indeterminate factor could 
be considered to be negligible. Otherwise, the contribution from the indeterminate factor could be significant. To examine 
this issue, we defined an angle (^9 ) at the z th iteration between search directions generated using determinate ^determinate 

and analytical sensitivities as 


0 Z - = cos 

The angle would be zero (Q f = 0) if the determinate and the analytical gradients were identical; otherwise, it would be 

nonzero ^ 0) . In some scale, the angle was a measure of the difference between the closed-form and determinate 

gradients of the active constraints. For a five-bar truss, the angle 0/ was calculated at the initial and the optimum design 
points for the SQP algorithm. Numerical values for the directions and angle are given in table 14. 


^anl ' d determinate 


*anl 


d, 


(34) 


determinate 
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Table 14. Angle between search directions solved with SQP algorithm 

for the five-bar truss 


Direction vector, d 

Initial design point, 

Optimum solution point, 

angle ©initial = 25° 

angle 0 op t = 64° 

Sensitivity 

Sensitivity 

Determinate Analytical 

Determinate Analytical 


initial 

a determinate 


' 0.38 


0.31 


1.51 



-1.49 

-0.16 


-0.03 


-1.44 



0.19 

-0.23 

a*. 

B.B. 

p] 

II 

-0.36 

> d opt 

determinate 

4.14 

■KT 5 

{rf} anl -° pt = . 

1.77 

-0.39 


-0.22 


-0.30 



-0.30 

-0.09 


-0.07 


0.00 



0.00 


The angle was different at the initial as well as at the optimal design points with values ^^1 = 25° and op t = 64° . For the 
three large problems, the angles were calculated for the SQP algorithm. At the initial point, the angles varied from 20° for 
the forward-swept wing to 51° for the 20-bay truss. At the optimum solution point, the minimum and maximum values for 
angles were 39° and 89°, respectively. In summary, the directions taken from the initial point to reach the optimum solution 
were different for the analytical and the determinate sensitivities. In other words, the indeterminate component of the 
sensitivities in was not negligible and it changed the path of optimization. 

Design optimization was not sensitive to the precision of the gradients of the stress and displacement constraints of an 
indeterminate truss. Simple determinate sensitivities performed very well for all design problems. The CPU time to optimum 
solution was substantially reduced with the determinate sensitivity. The question is: Why did the determinate sensitivity 
outperform the full closed-form gradient? An answer was attempted through a discussion of the nature of the structural 
design problem. 


4.4 Nature of the structural design problem 

The nature of the design optimization problem was examined by considering the three-bar truss. The problem had three 
design variables (A h A 2 , A 3 ). It had three stresses ( G x , <J 2 , <T 3 ) and three stress constraints. Likewise, there were two 
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displacements (X h X 2 ) and two constraints for each load case. There were three implicit structural analysis relationships 
between the five behavior variables: 


0 "l — C >2 "T CJ 3 — 0 

(35a) 

X l =^t ( ct 2 - 2 CT l) 

(35b) 

X 2 = —{ a 2) 

(35c) 


Equation (35a) is the compatibility condition and it is expressed in stresses. Equations (35b) and (35c) are the deformation 
displacement relations, also expressed in stress variables. It is important to observe that the three implicit relationships do not 
explicitly contain the design variable. In other words, gradient of one stress can be expressed in terms of the gradient of other 

f £ \ 

stresses. Such as for example, (V<T 2 = V [<J X + <r 3 )) and VX, = — (V<T 2 — 2 V <T, ) .If the second stress and second 

displacement constraint become active, then the coefficient matrix [Q\ of the direction vector {d} will become singular. 
Singularity will be avoided when determinate sensitivity is used because it restores the full rank of the sensitivity matrix. In 
other words, the implicit relationship of the behavior constraints induced singularity in the design optimization of structures. 


4.5 Coefficients in Equilibrium and Compatibility Matrices 

The number of stress and displacement components in the implicit relationship, similar to that in equation (35), depends on 
two quantities: q ee and q cc . The number of entries (or nonzero coefficients) in a column of the equilibrium matrix [ B ] is q ee . 
Likewise, the number of nonzero coefficients in a row of the compatibility matrix [C] is q cc . Both q ee and q cc are small 
numbers. In other words, a small number of stresses are dependent. Likewise, a displacement is dependent on few stresses. 
Consider the example of the 20-bay truss. The number of coefficients in a column of its equilibrium matrix varies between 
one and four (1 < q ee < 4; four is more prevalent). The coefficients in a row of the compatibility matrix range between 6 and 
20 (6 < q cc < 20; the typical number is six). Six typical dependence relationships for the truss are discussed next. 


For case 1, one stress and one displacement are dependent because q ee = 1. For case 2, one stress is dependent on two 
displacements because q ee = 2. For case 3, one stress is dependent on four displacements because q ee = 4. For cases 4 and 5, 
six stresses are dependent because q cc = 6 . For case 6 , 20 stresses are dependent because q cc = 20. The sixth case is interesting 
because the 20 stresses belong to the bottom chord members of the truss. The bottom chord is a natural load path but it can 
promote singularity in the optimization process, which however, can be avoided when approximate sensitivity is used. For a 
truss, a small number (in the range of two to six) of active stress and displacement constraints can be dependent. 


A traditional structural optimization problem contains dependent constraints. A small number of active stress and 
displacement constraints can be dependent. The multitude of implicit constraints reduces the rank of the coefficient matrix 
\Q]. Simple determinate sensitivity worked well because it restored the full rank. We suggested (ref. 22) that a set of 
independent constraints should be separated out of the given stress and displacement constraints by using a singular value 
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decomposition algorithm. The exercise has to be performed before the generation of each search direction. This technique 
works well for small problems. For larger problems, the decomposition process increases the numerical burden in 
optimization, which is already computationally intensive. The current recommendation is to use simple determinate 
sensitivity because it restores the full rank of the matrix [ Q ] and converges more rapidly. 


4.6 Adjustment for the stiffness method 

Design is stress driven both for stress and displacement limitations. The area of a truss bar for stress limitation can be updated 


as 


A n 


= A old . Likewise, the displacement formula (X = JGF ) can be manipulated to obtain an area update formula 


for the stiffness limitation (ref. 26). The two features make the method of force an attractive tool for design applications. We, 
however, realize that the stiffness method is very popular. The question is, “Can sensitivities be approximated when the 
stiffness method is used as the analysis tool in optimization?” Such an approximation is straightforward for the stress 
constraints. It may pose a challenge for the displacement limitation because it is a global variable. In the stiffness method, the 
stress sensitivity can be obtained by dividing the force or stress parameters by the square of area or area, respectively, ( -F/A 2 ), 
The force or stress output of a stiffness code can be adjusted to obtain the approximate sensitivity for the stress constraints. 


The difficulty encountered in calculating the approximate displacement sensitivity is illustrated by considering the three- 
bar truss as an example. Consider one term in the derivative of the first displacement X x with respect to the area A\, For 
simplicity, one load component is set to zero, P y = 0. The closed-form derivative can be written as 


dX 1 

~dA x 


+ 4A 2 A 3 + 2^2A^£P X 

2 

e(^ 2A 2 A 3 + 2A } A 3 + 4lA x A 2 ) 


(36) 


The procedure of separating the sensitivity into determinate and indeterminate cannot be directly extended to the stiffness 
method. The stiffness method, in general, appears to have two major impediments for design calculations. 

(1) The method has fewer equations m than the number of design variables n: m < n. The three-bar truss has three 
design variables, three bar stresses, but there are only two stiffness equations. Three equations are required to size 



F i 


F? 

, and A 3 = 

f 3 


the three bar areas, like A± = 

,a 2 - 

Z 

J 

. In other words, it is not easy to link bar areas to 

CT 0 


CT 0 




displacements. At best, this link would provide a relationship of three design variables to two displacements, and it 
would not be a one-to-one mapping. 


(2) The stiffness method does not allow free movement between analysis variables like IFM, which 
allows movement from force to displacement { X} = \J][G]{F}, and vice versa, {F} = [G] -1 [#] r {A}. The 
two formulas, along with their governing equation, [S]{F} = {P}, make IFM a very attractive tool for 
sensitivity calculation and design optimization. 
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4.7 Extension to other structure types 

Extension of the approximate expressions is straightforward for beams and framework. Consider, for example, a beam with 
moment as the analysis variable and moment of inertia / and depth d as the design variables. The flexure formula can be 
differentiated to obtain the approximate sensitivity for stress constraints: 


<J = 



Sa- 

'di 


Md 


2 1 2 


- + 


neglect — » 


d dM\ 

y21~df, 


( 37 ) 


The derivative of the moment in the displacement formula (X = JGF ) would provide the approximate sensitivity of the 
stiffness constraints. In other words, in IFM, the generation of the closed-form and approximate expressions for beams is 
quite straightforward. The logic can be extended for framework that would require a mixing of the truss and beam 
expressions. Then computer software has to be developed to compare design optimizations using approximate and analytical 
sensitivities for such structures. The exercise is worth the effort because singularity can be eliminated to make optimization 
robust for flexural structures. 


4.8 Generalization of sensitivity approximation 

It is straightforward to extend the sensitivity- approximation concept to a general type of structure. For illustration, consider 
the example of a stress constraint, g(x) which is a function of design variables (x). The constraint can be expressed as a 
product of an explicit function F(x) and an implicit function R(x) as:. 


gix) = F (x).R(x) 


Its gradient can be written as: 


Vg(x) 


simple-Retain 


RVF 


calculation-mi ensive-DROP 

FVR 


simple-Retain 


RVF 


( 38 ) 


( si mple - Re ta in \ 

RVF ) in the gradient expression. This term is simple to calculate. Use of this 
simple term eliminated singularity in structural optimization because it has a full row rank. The recommendation is to drop 


the second term 


calculation - int ensive-DROP 

FVR 


because its generation is computation intensive and its retention reduces the rank of 
the sensitivity matrix, which makes structural optimization a singular problem. Consider a plate flexure problem as an 
example. Its thickness (h) is the design variables and it has two principal moments and M 2 ) . The von Mises stress, 

[a j ) that can be used to define the stress constraint, can be written as: 
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The gradient of von Mises stress can be obtained as: 


simple-retain 



The gradient is simplified by retaining the first simple term which can be generated with a trivial amount of computation 
as: 


simple-retain 



The approximate gradient concept is independent of structure type and it can be used in finite element analysis. Likewise, the 
gradient of displacement can be approximated for a general type of structure through the flexibility matrix. 

There are numerous dependent constraints in the design optimization problem of an indeterminate truss. A small set of 
stresses can be dependent. A stress also can depend on a few displacements. A truss design with many sets of dependent 
constraints may not be a well-posed mathematical programming optimization problem. However, it is a real-life industrial 
design problem. In optimization calculations, all constraints should be used in defining the feasible region. Sensitivities of 
only the independent constraints should be used to calculate the direction vector. The independence criterion will be satisfied 
when the proposed simple determinate design sensitivity is used. The optimum solution was reached with the determinate 
sensitivity even though the search directions differed for the determinate and analytical sensitivities. The use of simple 
determinate sensitivity substantially reduced the CPU time to solution. The integrated force method is an efficient analysis 
tool for the calculation of determinate sensitivity in particular and for design application in general. The concept of using 
approximate sensitivity in design optimization should be extended to flexural structures like beams and framework. 

5.0 Summary 

The cascade optimization strategy solved the subsonic aircraft design optimization problem, even though restarts 
were required. It is preferable to restrict the behavior parameters from approaching zero values. The optimum 
aircraft weight calculated by the Flight Optimization System (FLOPS) analyzer and the regression method 
approximation matched well. The deviation in the design variables between the two analyzers was not significant. 

Deviation can be significant for some behavior constraints when these constraints approach zero values. Overall, the 
performances of the neural network and regression method were comparable. The neural network followed a mean 
path, while the regression method exhibited a tendency to closely follow the FLOPS solution. For a single analysis 
cycle the FLOPS time measured in seconds was reduced to milliseconds by the approximators. The training, 
validation, and solution required a small fraction of FLOPS analysis and design time. For design optimization, the 
central processing unit (CPU) time with the FLOPS analyzer measured in hours was reduced to minutes by the 
neural network, and seconds by the regression method. Generation of high-fidelity input/output pairs to train the 
approximators was time consuming. 
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Chapter 5 


Stochastic Analyses and Optimization 


1.0 INTRODUCTION 

Traditional analysis and design of structures assumed parameters like the load, material 
properties, sizing variables and the like to be well defined or deterministic in nature. For example a steel 
column with a circular cross section with a six inch radius can be analyzed for 10 kip load, with a 
Young’s modulus of 30 ksi and strength of 20 ksi. However, engineers have recognized the existence 
of scatter or uncertainty in the material properties, in load, and in design parameters. Consider for 
example the yield strength of steel that is required to design a steel structure. Strength is measured in the 
laboratory from tests conducted on standard coupons. It is commonly observed that repeated tests yield 
different values for the strength of steel. The test data can be processed to obtain a nominal or mean 
value and a dispersion range or a standard deviation. The nominal strength along with a safety factor is 
used to define allowable strength for traditional deterministic design calculations. Alternatively, strength 
can be considered as a random variable with a mean value and a standard deviation. The experimental 
data can be processed to obtain a probability distribution function such as, for example, the commonly 
used normal distribution function that is defined by a mean value and a standard deviation. This concept 
for strength can be extended to Young’s modulus, Poisson’s ratio, density, and so forth and a probabilistic 
material model can be generated. The procedure can be repeated and a probabilistic load model can be 
developed for mechanical, thermal and initial deformation loads. Likewise, a probabilistic design model 
can be developed for sizing variables like depth and thickness of a beam. 

Because of the stochastic nature of load, material properties, and sizing variables, structural response 
consisting of stress, strain, displacement, and frequency become random parameters with mean values 
and standard deviations. The cumulative distribution concept can be utilized to estimate the value of a 
response parameter for a specified level of probability. For example, the value of displacement in a 
particular structure at a location in a direction can be estimated to be less than or equal to 0.18 in. for a 
25-percent probability of success or reliability. The value can change to 0.23 in. for 75 percent 
reliability. The concept illustrated for displacement can be extended to other failure modes or design 
constraints, which become a function of reliability. In other words, a structure can be designed for a 
specified reliability between 0 and 1 . High reliability can lead to a heavier design. The design is likely to 
be lighter when reliability is compromised. In other words, the weight of a structure becomes a function 
of the reliability. It will be shown that reliability versus weight traced out an inverted-S-shaped graph. 
Reliability can be changed for different components of a structure, like that of an airframe. For example, 
the landing gear can be designed for a very high reliability, whereas it can be reduced to a small extent 
for a raked wingtip. 
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A stochastic design optimization methodology (SDO) has been developed to design components of an 
airframe structure, like the Boeing 767-400ER extended range airliner made of metallic and composite 
materials. The design was obtained as a function of the risk level, or reliability p. The design method 
treated uncertainties in load, strength, and material properties as well as design parameters as 
distribution functions, which were defined with mean values and standard deviations. A design 
constraint or a failure mode was specified as a function of reliability p. The SDO capability has been 
obtained by combining three codes: 

1. The MSC/Nastran code (Ref. 1) was the deterministic analysis tool. 

2. The fast probabilistic integrator, or the FPI module of the NESSUS software (Ref. 2), was the 
probabilistic calculator, and 

3. NASA Glenn Research Center’s optimization testbed CometBoards (Ref. 3) became the 
optimizer. 

The SDO capability requires a finite element structural model, which can be deterministic in nature 
because the coordinates of the grid points are not altered. Probabilistic models are required for material 
properties, the load cases, and design variables. Consider area of a bar with a circular cross section with 
a radius R. Its probabilistic model considers the radius as a random variable with a mean value (// s ) 

and a standard deviation (<r fl ) . The standard deviation expressed as a percent of the radius is a 

manufacturing issue and it need not be changed during theoretical design calculations. The stochastic 
optimization concept is illustrated considering an academic examples and a real-life raked wingtip 
structure of the Boeing 767-400 extended range airliner made of metallic and composite materials. 

Reliability-based design optimization requires a probabilistic analysis tool. Several such available tools 
are discussed in references 4 and 5. The quadratic perturbation technique should be the most preferred 
probabilistic analysis method because it provides closed form formulae for the mean value and standard 
deviation for response parameters like, stress, displacement and frequency. At this time, the formulae, 
have been programmed, for academic problems. For industrial strength problems, like the airframe 
wingtip structure, the Fast Probabilistic Integrator (FPI) module of the NESSUS code was used 
for probabilistic calculation. The probabilistic response is used to formulate the stochastic 
design problem, ft is solved using the optimization test-bed, CometBoards. The probabilistic 
analysis and design concepts are illustrated for an academic example and for an industrial 
strength aircraft wingtip structure made of metallic and composite materials. 


NASA/CR— 201 3-2 1 7748 


164 



2.0 Probabilistic Structural Analysis 

Popular probabilistic analysis formulations include Monte Carlo simulation (ref. 6), sampling 
and stratified sampling techniques (ref. 7), the Latin hypercube technique (ref. 8), response surface 
method, and others. Monte Carlo simulation is a powerful numerical approach, but it is repetitive and 
computationally expensive. Numerical integration, second-moment analysis, and stochastic finite 
element methods are also available. The perturbation method has been used extensively in developing 
the stochastic finite element method because of its simplicity, efficiency, and versatility. Cambou (ref. 9) 
proposed a first-order perturbation method for the finite element solution of linear static problems. 
Cornell (ref. 10) introduced a second-moment concept. During the 1980s, the method was developed 
further by Hisada and Nakagori (ref. 1 1) for static and eigenvalue problems. The perturbation method 
was also adopted by Handada and Anderson (ref. 12) for static problems of beam and frame structures. 
Applications of perturbation methods have also been reported by Benaroya (ref. 13), Eleishakoff (ref. 
14), and Schueller (ref. 15). For academic problems a quadratic perturbation technique is employed to 
calculate the mean value and the covariance matrix in closed form for stress and displacement. 

Stochastic analysis formulas have been programmed analytically in Maple V software, also in 
FORTRAN language. Solutions for simple structure have been obtained and verified 
via Monte Carlo simulations. For the wingtip industrial problem, fast probability integration (FP1) is 
used, which is based on a most probable point (MPP) concept. The majority of the probability is 
concentrated near the MPP. The MPP is located and different probabilistic methods are used to 
determine the probability in the failure prone region. The probabilistic response is generated from a few 
deterministic solutions. 

2.1 Stochastic Sensitivity 

The sensitivity analysis of structural systems with respect to variations in adjustable or design 
parameters is one of the approaches to evaluate the performance of structures. The process identifies the 
set of critical design parameters. In general, sensitivity analysis of structural system involves 
computation of the partial derivatives of some response function with respect to the design variables. 
Parameter identification and reliability assessment is important for system optimization. 

However, the conventional sensitivity analysis of structures is based on the assumptions of complete 
determinacy of structural parameter. In reality the occurrence of uncertainty associated with the system 
parameters is inevitable. Thus, it is necessary to estimate the effect of uncertainty in stochastic 
sensitivity derivatives with respect to random design variables. The stochastic structural sensitivity 
analysis is concerned with the change of stochastic structural response due to variations in stochastically 
for described design parameters. There are several methods for calculating the sensitivity of structural 
response to change in design parameters, including analytical, numerical and semi-analytical methods. 
Semi-analytical method is employed in the stochastic sensitivity analysis for Integrated Force 
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Method/Dual Integrated Force Method (IFM/IFMD), since this method based on the perturbation 
technique is easy to implement as a numerical approach and it is an efficient analytical approach. 


2.2 Perturbation Method 

The perturbation method for probabilistic analysis has been developed for both force and 
displacement methods. The method of force used is referred to as the integrated force method (IFM). 

The dual of the primal IFM, or the dual integrated force method (IFMD) became the displacement 
method. The perturbation technique yields the mean value, and the covariance matrix for the response 
variables like internal force and displacement. The cumulative distribution concept is applied to 
determine the response for a specified distribution function and a probability level. The parameters for 
perturbation, or the primitive random variables, are separated into three groups. These are loads, 
material properties and design variables. The load vector encompasses mechanical load, thermal load, 
and settling of support or initial deformation load. It is defined in terms of a mean value for the load 
vector and a covariance matrix. Material properties considered are elastic modulus, Poisson’s ratio, 
coefficient of thermal expansion, and density. Their mean values and the associated covariance matrices 
are the stochastic parameters. The sizing design parameters are the areas of bars, moment of inertia of 
beams, thickness of plates, and so forth. Their stochastic parameters are defined through mean values 
and a covariance matrix. The geometrical parameters like length of a truss bar and spans of a plate are 
considered as deterministic variables. 

Dr. Wei’s doctoral thesis (ref. 16) is devoted to stochastic analysis and design optimization for 
structural system. In this thesis, a full set of stochastic analysis formulas for IFM/IFMD has been 
derived by using the first- and second-order moment perturbation methods. The sensitivity analysis of 
stochastic responses has also been studied. Neumann expansion technique for stochastic analysis using 
IFM/IFMD has also been discussed. Furthermore, the variational energy formulation for stochastic 
analysis for IFM is provided. Those closed-form expressions for stochastic analysis have been 
programmed in Maple V and a FORTRAN code was also developed. Solutions were obtained for a set 
of one dozen examples. The closed form solutions have been verified and compared with Monte Carlo 
simulation, Neumann expansion technique using ANSYS software (ref. 17). Stochastic optimization 
was formulated for stress and displacement constraints for a specified probability of occurrence. The 
objective function was the weight of the structure which was also a function of probability. The mean 
value and/or standard deviation were considered as design variables. A set of more than one dozed 
problems were solved. Solution of weight versus reliability produced the inverted S-shaped graph for the 
problems. As mentioned earlier, the center of the inverted-S graph corresponded to 50 percent (p = 0.5) 
probability of success. A heavy design with weight approaching infinity could be produced for a near- 
zero rate of failure that corresponds to unity for reliability p(orp = 1). Weight can be reduced to a small 
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value for the most failure-prone design with a reliability that approaches zero (p = 0). Reliability can be 
changed for different components of an airframe structure. For example, the landing gear of an airliner 
can be designed for a very high reliability, whereas it can be reduced to a small extent for a raked 
wingtip. The SDO capability is obtained by combining three codes: The MSC/Nastran code was the 
deterministic analysis tool, The fast probabilistic integrator, or the FPI module of the NESSUS software, 
was the probabilistic calculator, and NASA Glenn Research Center’s optimization testbed CometBoards 
became the optimizer. The SDO capability requires a finite element structural model, a material model, a 
load model, and a design model. The stochastic optimization concept is illustrated considering an 
academic example and a real-life raked wingtip structure of the Boeing 767-400 extended range airliner 
made of metallic and composite materials. 

In the perturbation technique, the mean value and the covariance matrix of a response variable are 
obtained in two basic steps. First, a response variable is expanded in a Taylor’s series with respect to the 
primitive random variables. The linear and quadratic terms are retained in the series expansion. An 
application of the expectation operator E[ ..] yielded the expressions for mean value and the covariance 
matrix. Prior to the expansion, the primitive parameters are scaled to obtain a normalized variable with a 
zero mean and a small variance. In other words a two random variable problem is transformed to a 
single normalized random parameter with a zero mean value. 

Consider, for example, the coefficient of thermal expansion a with a mean value of ju a and a standard 
deviation of G a . Its normalized variable q a is defined as 


9 a = 


a -Pa 
Va 


( 1 ) 


The mean value of q a is obtained as zero by using the expectation operator (E[..]): 


Na E 


a-Va 

Va 


E[«]-E [n a ] 




f^a L L a 


= 0 


(2a) 


The variance of q a is obtained as 


<x 2 ^ 

A 


(2b) 


The variance of a is assumed to be much smaller than the square of its mean. This justifies the use of a 
Taylor expansion of response variables in terms of the normalized primitive random variables. The 
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Taylor series expansion is discussed considering the example of the internal force vector {F}. The 
expansion can be expressed as 


{F} = {F} + ({^} r £{F}) + i^} r C{F}C T {?}) + h.o.t 


(3) 


where 


£ (-))-E 


<li 


%• 




and 


{{q} T £(•)£ 

V J Uj 


li ± 

dqjdqj 




Similar expansions are also obtained for the governing matrix [S] of IFM, the load vector {P}, 
displacement {X}, the flexibility matrix [G], and initial deformation {/? }, where their deterministic (or 
over-bar) quantities are calculated by setting the normalized primitive random variables to zero (q t = 0). 
Four steps are followed to derive the expressions for the mean value and the covariance matrix of the 
force vector: 

Step 1. Substitute the Taylor series expansions into the IFM governing equation. 

The left-hand side contains a single zero-order term, two linear terms, and three quadratic 
terms in the normalized primitive random variables, and the right-hand side has one term of 
each order. 

Step 2. Equate equal order terms. The zero-order term gives the deterministic solution. The linear 
terms drop out. The second-order terms produce an equation with four quadratic terms on the 
right-hand side and a single term on the left-hand side. 

Step 3. Take the expectation E[..] of terms in Step 2 to obtain the mean value of the force vector. 

Step 4. Take the expectation E[..] of the product {F\ {F} T , and subtract the square of the mean of 

the force vector, as calculated in Step 3, to obtain the covariance matrix of the force vector. 

The expressions for mean value and covariance matrix for force and displacement vectors via 
integrated force method, is not reported here. 
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The derivation included mechanical and initial deformation loads. The mean value of the force and displacement 
contained zero-order and quadratic terms, but not linear terms, in the Taylor series. However, the 
covariance matrix retained the zero-order, linear, and quadratic terms of expansions. For simple 
academic structures the expressions were programmed in closed form in Macsyma software (ref. 18). 


3.0 Literature on Stochastic Optimization of Structures 

In the past decades, optimization under uncertainty has progress rapidly, both in the theoretical 
formulation as well as in the development of algorithms. The approaches included: expectation 
minimization; minimization of deviations from goals; minimization of maximum costs and optimization 
over soft constraints. Introducing uncertainty within the deterministic optimization framework allows for 
the more realistic results, but also adds the complexity to optimization. There are two fundamental 
techniques to solve the stochastic nonlinear optimization problem; The first one is the regularized 
decomposition. Ruszczynski (ref. 19) added a quadratic terms to the objective function for improving 
convergence of the L-shaped decomposition method. Taking advantage of the augmented Lagrangian 
method, Dempster (ref. 20) added a quadratic penalty to ensure convexity, yielding more efficient 
computation. Rockafellor and Wets (ref. 21) also developed a similar method, called the Progressive 
Hedging Algorithm. Their significant limitations are on either the objective function type or the 
underlying distributions for the uncertain variables. The other technique is to identify problem for 
specific structures and transform the problem into a set of deterministic nonlinear programming 
problem. For example, Chames and Cooper (ref. 22) replaced the uncertain constraints with the 
appropriate probabilities expressed in terms of moments. In chance constrained programming, the 
uncertainty distributions should be stable; the uncertain variables are linear in the chance constraint, and 
the problem need to satisfy the general convexity conditions. Recently Marti (ref. 42) converts the 
random problems in plastic structural analysis and optimal design into deterministic substitute problems 
by applying stochastic optimization methods. The advantage of these methods is that one can apply the 
deterministic optimization techniques to solve the stochastic optimization problem. In the deterministic 
nonlinear optimization, the sequential quadratic programming (SQP) algorithms are widely considered 
today as the most effective general techniques for solving nonlinear programming with nonlinear 
constraints numerical experiments also show that SQP methods have good convergence properties. The 
idea of solving nonlinear problem by a sequence of quadratic programming subproblems was first 
suggested Wilson (ref. 23), SQP methods were popularized mainly by Biggs (ref. 24), Han (ref. 25) and Powell 
(ref. 26 and 27). Not only is SQP very efficient for solving nonlinearly constrained optimization 
problems, but also has been applied to some problems such as non-smooth equations, variational 
inequality problems, mathematical programs with equilibrium constraints, etc. and some large-scale 
problem. 


NASA/CR— 201 3-2 1 7748 


169 



Since the stochastic optimization problem is essentially a nonlinear program, the solution methods for 
nonlinear programming problems can be applied, and some recent studies indicated that sequential 
quadratic programming (SQP) is a very promising method particularly for solving nonlinear programs. 
Moreover, the SQP method had been successfully adapted for solving optimal design of structures based 
on the perturbation method. Lee and Lim (ref. 28) presented a general formulation of the optimism 
design problem with the random parameters. Sedaghati and Esmailzadch (29) developed a new 
structural analysis and optimization algorithm to determine the minimum-weight of structures with the 
truss and beam-type, members under displacement and stress constraints, by using the force method. 
The algorithm is based on the sequential quadratic programming (SQP) technique, and the finite element 
technique is base on the integrated force method. Chen (ref. 30) proposed a stochastic optimization 
modeling procedure to improve product quality by reducing variations. He developed two stochastic 
versions of SQP respectively embedded with a Monte Carlo simulation and numerical approximation in 
dynamic-characteristic robust design. Fares, et. al. (ref. 31) provided a sequential semi-definite 
programming, an extension of SQP, to solve the robust control problems. To improve the quality of a 
product caused by variation without eliminating these causes, the robust design is introduced to reduce 
the effects of variability. Doltsinis and Kang (ref.32) presents a formulation for the robust design of 
structures with stochastic parameters. Both the expected value and the standard deviation of the 
objective function are to be minimized in the optimization. The robust design optimization problem has 
been solved with SQP. Lee and Park (ref. 33) provided a multi-objective function defined by the mean 
and the standard deviation of the original objective function. To obtain the optimum value insensitive to 
variations of design variable, the constraints are supplemented by adding a penalty term to the original 
constraints. Levi et. Al (ref. 34) introduced a multi-objective optimization method suitable for dealing 
with stochastic systems. The optimization method can manage the uncertain system parameters and 
external disturbances introduced by means of a Gaussian stochastic process. Parkinson et. al. (ref. 35) 
proposed a general approach for robust optimal design and addressed the design feasibility and the 
control the transmitted variation and also discussed the feasibility robustness and sensitivity robustness 
for robust mechanical design. Sundaresan et. al. (ref. 36) applied a sensitivity index optimization 
approach to determine a robust optimum. Mulve et. al. (ref. 37) treated robust design as a bi-objective 
non-linear programming problem and employed a multi-criteria optimization approach to generate a 
complete and dfficient solution set to support decision-making. Chen et. al. (ref. 38) developed a robust 
design methodology to minimize variations caused by the noise and control factors. The various 
objective functions, analysis techniques and applications used for robust design are reviewed by Zang 
(ref. 40). Mattson and Messac (ref. 41) also give a brief survey on how various robust design 
optimization methods handle constraint condition. 
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4.0 Stochastic Response Analysis 

Stochastic response analysis is illustrated considering the example of a three-bar truss shown in figure 1 
The geometrical dimensions of the structure like the coordinates of the four nodes and bar lengths are 
considered deterministic in nature. The random variables included mechanical load, thermal load, load 
due to settling of supports, Young’s modulus, coefficient of thermal expansion and bar areas. The truss 
has a total of 10 primitive random variables for the material properties, load, and sizing design 
parameters. 


4.1 Material Model 

The truss is made of steel, and it has two random material variables: the Young’s modulus and the 
coefficient of thermal expansion. The Young’s modulus E has a mean value of //£■ = 30 000 ksi and 
standard deviation of 3000 ksi. The coefficient of thermal expansion a has a mean value of //„ = 6.6x 10“ 
6 /°F and standard deviation of 6.6x10 7 /°F. 


100 in. 100 in. 
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4.2 Mechanical Load Model 


The structure is subjected to mechanical and thermal loads as well as loads due to settling of support. 
The entire structure is subjected to a single change of temperature T. Both mechanical and thermal 
loads have two components and covariance matrices are considered to define the loads. Mechanical load 


is defined by its mean value and a covariance matrix. The load 



is applied at the free node 1 . 


Mean value: 




kip 


Covariance matrix: cov 

V 


Px 

P Y 


25.00 6.25 

6.25 25.00 


4.3 Thermal Load Model 


Node 1 of the structure is subjected to a change of temperature, with a mean of 
ju T = 50 °F and a standard deviation of 5 °F. The other three nodes are at ambient temperature. Thermal 
load is obtained for the bar temperature, which is calculated from the average of the temperatures at its 
two connecting nodes. 


4.4 Load due to settling of Support 

Support node 2 settles in both the x- and y-directions by amount J ^ 2X 1 inch as shown in Fig. 1 . 

l A 27 J 

Settling is accommodated by a mean value and a covariance matrix. 


Mean value: 



= in. 


Covariance matrix: 


cov 


x 2X 
^2 Y 


25.0 37.5 

37.5 225.0 


x!0“ 
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Covariance matrix: cov 
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The 10 random variables of the problem are scaled to obtain the 10 primitive variables ql, q2,...,ql0, 
which are defined as 


4 
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The normalized primitive random variable, {q}, with zero mean, and standard deviation given by the 
ratio of the standard deviation to the mean of the corresponding stochastic variable is assumed to be of 
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order o(l). This justifies the use of a Taylor series expansion in {q}. Only up to quadratic terms are 
retained in the series. 


The deterministic response for forces and displacements are as follows: 

' 62.76' 

Bar forces: {f} = {^i f } = \(>\.2 a\ kip 

- 7.95 


Nodal displacements: {x} = {n x ) 


0.201 1 
fin- 

-0.24lJ 


The probabilistic response for bar forces consists of the mean values and the covariance matrix as 
follows: 


Mean values of bar forces: { uy } = 


62.78 

61.22 

-7.93 


kip 


Covariance matrix of forces: cov(f) = 


19.27 

1.78 

-4.48 

1.78 

18.93 

-8.83 

-4.48 

-8.83 

21.78 



2 


1.0 

0.093 

-0.219 

Correlation matrix for forces: 

pa 

- PF = 

0.093 

1.0 

-0.435 


_( pi\pj)_ 

-0.219 

-0.435 

1.0 
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Likewise, the probabilistic response for nodal displacements are 


Mean values of displacements: {^} 


0.197 

- 0.237 


in. 


Covariance matrix of displacements: cov(v) 


12.2 

- 8.9 


- 8.9 

8.7 


xKT 4 


Correlation matrix for displacements: Px = 


1.0 

- 0.866 


- 0.866 

1.0 


Where, the elements of Pf and p x are the correlation coefficients. 

For the bar forces, there is very little difference between the mean value { P p} and the deterministic 
solution )//,,} . For the nodal displacements also, there was little difference between the mean value {//j j and 
the deterministic solution { Px } . For all examples, big or small, it was observed that the mean value was quite 
close to the deterministic solution. 

The standard deviations for the bar forces were about 7 percent for the bars 1 and 2. However, it was 
59 percent for the bar 3. The variation in the primitive variables like mechanical, thermal, and settling of 
support loads was between 5 and 10 percent. The variation in the response or bar forces of the basic 
structure is comparable to that for the load variation. However, the variation in the redundant bar force 
at, 59 percent, is 6 times that of the load. 
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The standard deviations for the displacements are between 12 and 17 percent, which are higher than the load 
variation. The variations for the displacements are more pronounced than that for the bar forces because these are 
susceptible to loads as well as to material properties and design variables. The probability density function and the 
cumulative distribution function for the bar 1 force and magnitude of the vertical load are depicted in Figs. 2 and 
3, respectively. There is 50-perent likelihood that the first bar force (with a mean of 62.67 kips) is between 59.80 
and 65.73 kips. 

The value of the response variables for percent probability of success p = 50, 25, and 75 percent are listed in 
Table 1. A 5-percent change was seen in the forces of bars 1 and 2 between the p = 50 percent and p = 15 percent 
levels. For the bar 3 force, the changed noted was 40 percent. A 12-percent change was seen in the horizontal 
displacement between the p = 50 percent and 




Figure 3. Cumulative distribution function for bar 1 force and vertical load. 
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p = 75 percent levels. It was only 8 percent for the vertical displacement. 


Table 1. Stochastic response for three-bar truss. 


Response variable 

Upper limit of response variable (range) for probability of occurrence,/? 


50 percent (mean) 

25 percent 

75 percent 

Bar force, kip 




F i 

62.76 

59.80 (95%) 

65.73 (105%) 

F 2 

61.24 

58.30 (95%) 

64.17(105%) 

F 3 

-7.95 

-4.80 (60%) 

-11.09(140%) 

Displacement, in. 




u 

0.201 

0.178 (88%) 

0.225 (112%) 

V 

-0.241 

-0.221 (92%) 

-0.261 (108%) 
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Figure 4. Sensitivities for bar 1 force with respect to 10 random variables. 
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Figure 5. Sensitivities for vertical displacement with respect to 10 random variables. 


The sensitivity of the bar 1 force and the vertical displacement at the 75-percent probability of 
occurrence level with respect to the 10 primitive random variables are shown in Figs. 4 and 5, respectively. 
The bar force is equally sensitive to the areas of bars 1 and 2 in magnitude, but opposite in directions. The 
sensitivity of the bar 1 force with respect to the bar 3 area is very small. The vertical displacement is 
equally sensitive to the areas of bars 1 and 2, but is insensitive to the bar 3 area. The Young’s modulus 
has little effect on the bar force but has the most effect of all the primitive variables on the vertical 
displacement. The sensitivities with respect to thermal coefficient and temperature are identical because 
the analysis uses their product as a single entity. Both load components strongly effect the bar force, 
whereas only the vertical load effects the vertical displacement to the same degree. Temperature and 
settling of support loads have a moderate effect. 

5.0 Comparison With Other Methods 


Stochastic response via perturbation method was compared with other probabilistic formulations 
considering the example of the three-bar truss shown in Fig. 1 with the same 10 random variables for the 
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Material properties, load, and sizing design parameters. A summary of response calculated by the four Methods is 
depicted in Table 2. Two different computers were used in the calculations shown in Table 2. A Dell Inspiron 
desktop with four CPU and 3.2 GHz processor was used to obtain most of the results. The FPI code shown in the 
last column was run in a very fast Red Hat Linux dual CPU workstation with a 32 GHz processor. The four 
methods compared were 

(1) Perturbation method (PM) 

(2) Direct Monte Carlo simulation (DMCS) 

(3) Latin hypercube simulation (LHS) 

(4) Fast probabilistic integration of the NESSUS code (FPI) 


Table 2. Probabilistic response compared for three-bar truss. 


Parameter 

Perturbation method 
(PM) 

Direct Monte Carlo 
simulation (DMCS) 

Latin hypercube 
simulation 
(LHS) 

Fast probability 
integrator 
(FPI) 

Mean 

value 

Standard 

deviation 

Mean 

value 

Standard 

deviation 

Mean 

value 

Standard 

deviation 

Mean 

value 

Standard 

deviation 

f Fj ] Kip 
Force :< F 2 f 

LJ 

■< 

'62.77' 

61.24 

-7.95 

> 

< 

4.39 

4.35 

4.67 

> 

< 

'62.77' 

61.27 

-7.93 

> 

< 

4.39^ 

4.35 

4.67^ 

> 

< 

'62.76 

61.25 

-7.96 

> 

< 

4.39 

4.35 

4.67 

> 

< 

'62.78' 

61.21 

-7.93 

< 

4.39 

4.35 

4.67 

> 

r "I Ksi 

°i 

Stress • i cr 2 r 

l°3 J 

c 

'63.29 

61.70 

—3.98 

> 

< 

6.71 

6.20 

2.34 


< 

63.16 

61.58 

-3.97^ 

► 

< 

5.89' 

5.42 

2.79 

> 

< 

'63.15 

61.57 

-3." 

> 

< 

5.77' 

5.38 

2.75 

> 

< 

'62.78' 

61.21 

—3.96 

< 

6.34 

6.03 

2.69 

> 

r " | in. 

Displacement : \ > 


' 0.20 
[-0.24, 



0.0041 

0.003J 

l 

0.20' 

[-0.24J 


< 

0.0041 

0.003J 

l 

0.20 1 

[-0.24J 


< 

0.0041 

0.003J 

l 

' 0.201 
[-0.24J 


0.0041 

0.003J 

Computation Time 

CPU seconds 

7.0 

4245 

390 

1 

Normalized time 

1.0 

606 

56 

— 
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The stiffness method as implemented in the ANSYS code was used in the DMCS and in the LHS. 
Response was calculated via the primal and the dual integrated force methods (IFM and IFMD, 
respectively) in the perturbation and fast probability integration techniques. Both IFM and IFMD yield 
identical solutions for deterministic as well as stochastic calculations even though the structure of 
equations differed, at least in appearance. Response for the three bar forces, stresses, and displacements 
by the four methods (PM, DCMS, LHS and FPI) and time to solution (CPU seconds) are depicted in the 
Table 2. For bar forces the mean values and standard deviations were in good agreement for all four 
methods: the perturbation method, direct Monte Carlo simulation, Latin hypercube simulation, and fast 
probability integrator. The displacements showed an almost perfect match across the four methods. 
There was a minor mismatch among the methods for bar stress. The maximum difference was about 0.8 
percent in the mean value of the stress, whereas it was about 19 percent for standard deviation for the 
redundant member. Direct Monte Carlo simulation required 12 500 samples for convergence, whereas 
1000 samples were sufficient for the Latin hypercube simulation. The time to calculate the response was 
very small, between 1 to 7 CPU seconds by the perturbation method as well as by the fast probability 
integrator. The calculation time increased many times for the Monte Carlo and Latin hypercube 
simulations. Monte Carlo simulation required about 4245 s, which corresponded to 606 times that 
required by the perturbation method. Latin hypercube method took 56 times as long. Overall, the 
performance was satisfactory for all the methods. The Ph.D. dissertation of Wei provides merits and 
limitations of different probabilistic solution methods. 

6.0 Probabilistic Design Optimization 

Probabilistic design optimization accounts for the uncertainties in the model, like those in design 
variables, in loads, in strength and in material properties. The design is obtained as a function of the risk 
level, or reliability, p. A design with a 50-percent rate of success (corresponding to one failure in two 
samples, or N = 2 and p = 0.5) is the mean-valued design. For the most reliable design the parameter p 
can approach unity. For example, for one failure in 50 000 samples,/? = 1 - 1 /N= 0.99998. For a design 
that is not reliable or prone to failure, p can approach 0, (e.g., for 999 failures in 1000 samples, p = 
0.001). Consider a deterministic design that can be obtained for the mean values of random parameters, 

like loads, materials, strength etc. Let this design be designated by |A det |. Consider next the mean 
value of the stochastic design { p x } . that will be discussed next. From the solution of more than two 
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dozen stochastic design optimization problems it was observed that the mean valued design was quite 
close to the deterministic solution, A det j = { // v } j. 

The constraints that represent failure modes are defined for a specified rate of failure. A probabilistic 
design optimization problem can be formulated with two merit functions: 

(1) to minimize the mean value of the weight 

(2) to minimize the variation or the standard deviation of the weight. 

Consider for example, the three-bar truss shown in Fig. 1. The area or diameter of a bar becomes the 
design variable. Probabilistic calculation can consider the determination of mean value and a standard 
deviation in the diameter. However, the standard deviation in the diameter is a manufacturing issue. A 
design process uses available members with specified diameter or cross-sectional areas. Probabilistic 
optimization can be separated into reliability-based design and robust design. Reliability-based design 
minimizes the mean value of the objective function. Robust design simultaneously minimizes the mean 
value as well as the standard deviation of the objective function. This paper emphasizes reliability-based 
optimization. However, for completeness, the robust design concept is discussed for the three-bar truss 
shown in Fig. 1 . 

The mean value and standard deviations can be defined as 


Mean value: //y = E [f(xj] = , * 2 , . . )dx\dx 2 . . .dxfc 


Standard deviation: ay = E 


2 2 

(f(x)-///) ( f Cv) - /// ) p(x h x2,...x n )clx\dx 2 ...dx k 


(5a) 

(5b) 


where 

p(x 1? x 2v . .x k ) is the joint probability density function 
f(x) is the merit (or weight) function 
Xi are the design variables 
E[..] is the expectation operator 


NAS A/CR— 20 13-217748 


181 



It is difficult to define the joint probability density function for an industrial-strength structural design 
problem. Normal distribution can be assumed for the random variables, which are further considered to 
be independent. The robust design optimization problem can be reformulated as follows: 


Minimize both mean value /// and standard deviation of f = [//y,ay] 
for constraints E [g, (x)] + p i o{g i (x)) < 0 
within prescribed bounds x L <x<x v 


where, 

/ is a vector of two objective functions 

/// is the mean value 

Of is the standard deviation 

Pi is a constraints (g , ) satisfactory index 


Pareto optimality criteria is a classical technique to solve a multi-objective problem. Because of 
difficulty, however, the problem is transformed to a single composite objective function and solved. A 
composite objective is obtained by a linear combination of mean value and standard deviation with a 
weighting factor a in the range 0 < a < 1 . 
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Table 3. Robust design for the three-bar truss. 


Design parameters 

Weight factor, a 

0 

0.25 

0.50 

0.75 

1.00 

Area: - 

^2 

M 

■ , in 2 

< 

1.1042 

0.7359 

0.1000 

> 

< 

1.1042' 

0.7359 

0.1000 

> 

< 

1.1042' 

0.7360 

0.100 

> 

4 

1.1042 

0.7360 

0.1000 

k J 

> 

< 

1.1042 

0.7361 

0.1000 

> 

Weight: { mean },lbf 
( variance] 

1 

'70.48701 
8.4500 J 

< 

'70.48631 
8.4499 J 

~l 

'70.48701 
8.4500 J 

~l 

'70.48891 
8.4502 J 

1 

'70.49081 
8.4504 J 


* * 

where y i w and cr w are the mean and standard deviation of weight, respectively, when solved 
independently. The robust design of the three-bar truss shown in Fig. 1 uses the data given for its 
analysis, along with a weight density that has a mean value of 0.289 lbf/in with a standard deviation of 
0.005 lbf/in . The mean value and the standard deviation for the strength are 20 ksi and 2 ksi, 
respectively. The optimum solution is given in Table 3 for different values of the parameter a for a 50- 
percent probability (p = 0.5). Mean value and standard deviation versus the parameter a are plotted in 
Fig. 6 with an exaggerated scale along the y-axis. 
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Figure 6. Robust design optimization for three-bar truss. 
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The mean value and standard deviation versus the parameter a traced out the typical graph as shown in 
Fig. 6. The robust design corresponding to the intersection point MS represents a simultaneous minimum 
for the mean value of weight as well as its standard deviation. The intersection occurred for a = 0.58. 
Weight has a mean value of 70.488 lbf and a standard deviation of 8.45 lbf. Critics may label the robust 
design optimization calculations as academic because the range of variation is very small for the mean 
value of the weight 70.4863 to 70.4908 lbf as well as for its standard deviation 8.4499 to 8.4504 lbf. 
Robust design optimization is discussed further in reference 4. Reliability-based optimization is referred 
to as stochastic design optimization and is discussed next. 

7.0 STOCHASTIC DESIGN OPTIMIZATION (SDO) 

A design, in stochastic optimization, is obtained as a function of the risk level, or reliability, p. A 
constraint or a failure mode is specified as a function of p. Optimum weight versus p traced an inverted- 
S-shaped graph. The center of the inverted-S graph corresponds to 50 percent,/? = 0.5, probability of 
success. A heavy design with weight approaching infinity can be produced for a near zero rate of failure 
that corresponds to unity for p, or p= 1. Weight can reduce to a small value for the most failure-prone 
design {p = 0). The stochastic design optimization (SDO) capability is obtained by combining three 
codes: (1) The MSC/Nastran code is the deterministic analysis tool, (2) the fast probabilistic integrator, 
or the FPI module of the NESSUS software, is the probabilistic calculator, and (3) NASA Glenn 
Research Center’s optimization testbed CometBoards became the optimizer. The SDO capability 
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requires four models: a finite element structural model, material model, load model, and design model. 
This report illustrates the design capability by considering the example of a Boeing 767— 400 ER raked 
wingtip. Stochastic design, also referred to as reliability based design optimization has been applied to a 
large-scale model of a vehicle side impact. 

7.1 FORMULATION OF RELIABILITY-BASED DESIGN 


The design testbed CometBoards has been extended into the stochastic domain. The objective is to 
determine the mean values of the design parameters that optimize the mean value of a merit function 
(such as weight), subject to a set of constraints. Variance of the weight can be back calculated. For 
stochastic optimization, the constraints gj are derived from random response variables within prescribed 
random upper and lower bounds g 1 - and gj , respectively, with a specified percent probability 

{p between 0 and 1), as: P( gj < gj < gj )>p. A stochastic constraint (e.g., a stress limitation x < To) can 
be expressed as 


Stochastic 

g(p) 


Standard deviation 


Mean value 


1 

/'r. 

- 


Ja 2 +a 2 

-1 

+ 

®*(nV n Tl ° 


1 {P 

_ '“no 



L Pm J 


(7) 


where 

p is the probability of occurrence 
ju T 1 is the mean value of stress 

is the allowable value of stress 

0 is the parameter phi critical 
cr ri is the standard deviation of stress 

is the standard deviation of stress allowable 
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The first term with superscript ’’Mean value” resembles a deterministic stress constraint. It is written in 
terms of the mean values of the response variable, here stress T\ has a mean value p T . Strength or 

allowable stress t\q has a mean value of // r|o . The second term with superscript ’’Standard deviation” 

accounts for the standard deviations <r Tj and a T in stress and strength t\ and rio, respectively. The 

inverse of the cumulative distribution function for the standard normal at probability level p is <!>*(/?). 
The stochastic constraint reduces to the deterministic form, when the variation is set to zero. 


Lim< 


a Tl 

a m 0 

ju r , „ — » allowable 
/ *10 


f 




Mr j 




-1 

V 

^10 

J 


a 2 + <r 2 

+ <£>*(/>) — — — ^<0 


P, 


T 10 


Deterministic constraint 


Mr* 




T 10 




( 8 ) 


The derivation of the constraint expression in Eq. (7) is straight forward. Consider a stress constraint 
r l ~ r ° ) with a probability p of occurrence of p: 


p(|r|-r 0 <°)> j p (9) 

Considering the case when x is greater than zero, the difference between the two random variables is set to a new random 
variable 9C 


9? = r - Tq 


( 10 ) 


The new random variable is normalized to obtain a standard normal random variable ® with mean 0 and variance 1 : 


® = 




°9t 


( 11 ) 


Thus 
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( 12 ) 


p(iR<0)>/? 


P 


f ,, ^ 
o<_M 

v y 


2 P 


The probability p(0 < x) is the definition of the cumulative distribution function for the standard normal F®. Thus 


r<t> 


Mm 

ij 


^ P 


(13) 


The minimum value of ® at which the probability level, p , is satisfied is obtained from the inverse of the cumulative 
distribution function for the standard normal. This value is labeled ®*. 


But 


So, 


(p ) = V (/>) 


MM = Mt~Mt 0 

2 _ 2 2 
erg; - cr r + £7 




< — 




0 

2 2 
OV +cr 
r 0 


(14) 


(15) 


(16) 


Finally, taking the case x < 0, normalizing with respect to the mean of the allowable and combining gives 


NAS A/CR— 20 13-217748 


187 



<0 


r- 
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Ja 2 + a 2 

Mr 

-1 

+ 

Mr, 

1 

> 


g<0 


( 17 ) 


As mentioned earlier, the constraint has two parts. The first part is similar to a deterministic constraint 
specified in terms of the mean value of the variable. The second part contains the contribution from stochastic: the 
inverse function and square of standard deviations as well as the mean value. For p = 0.50, <f>* is zero, and the 
constraint is similar to a deterministic constraint. Observe that ti is a response variable that is obtained from 
stochastic calculations and is not the same as the solution to the deterministic system. 

In the above derivation, the two cases x > 0 and x < 0 were treated separately. In doing so, an 
assumption was made that p T > <j t . This is reasonable when the constraint is active. For p = 0.25, ®* is - 
0.6745, and the constraint is relaxed. If it were an active constraint, the stochastic process becomes 
equivalent to increasing the allowable value. This would have a tendency to generate a lighter optimum 
design. For p = 0.75, ®* is 0.6745, the constraint is tightened. If it were an active constraint, the 
stochastic process becomes equivalent to decreasing the allowable value. This would have a tendency to 
generate a heavier optimum design. The stochastic design will have a tendency to produce a lighter 
design when p is less than 50 percent, but a heavier design otherwise. The formulation of reliability- 
based design optimization is quite similar to that for the deterministic problem. The mean value of the 
design variable is generated to minimize the mean value of weight subjected to the stochastic constraints 
defined in equation (7). 

8.0 NUMERICAL EXAMPLES FOR RELIABILITY-BASED DESIGN OPTIMIZATION 

Reliability-based optimization is illustrated considering two components of an airframe structure. The 
first component is a web of an airframe stabilizer, and the second component is a racked wingtip. 
Examples are proprietary components of the Boeing Company and provided to us as a courtesy to 
advance reliability-based optimization concept for industrial structures. Sufficient information will be 
provided to illustrate the concept. 

8.1 WEB PLATE OF AN AIR FRAME STABILIZER 

The inverted-S-shaped graph in reliability-based design optimization is illustrated considering the 
example of a web plate of a airframe stabilizer. The steel web plate, which was part of an aircraft 
stabilizer structure, was about 100 in. long, 12 in. wide, and 0.072 in. thick. It was idealized as a planar 
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structure in the x,y-plane. The finite element model had 80 CQUAD4 and 3 CTRIA3 elements and 103 
nodes. Load was applied in the x,y-plane. Load was obtained by projecting the reactive load from the 
stabilizer into the plane of the web. The maximum von Mises stress for the initial design was 3227 psi 
with a maximum displacement of 41.5 milli-inch. Load, material properties, and design parameters were 
considered as random variables. The mean values were taken equal to their deterministic values with a 
10-percent standard deviation. For example, the mean value of Young’s modulus was 30 000 ksi, (which 
was its deterministic value) with a 10-percent standard deviation of 3000 ksi. The stochastic response 
calculation used three different probabilistic analysis methods. All three methods used MSC/Nastran for 
deterministic calculation. The FPI module was also used by all three methods. The methods used were: 


Method 1: This method used the fast probabilistic integrator with normal distribution for all 
variables. 

Method 2: A neural network and a regression method were trained for the three design variables. 
Normal distribution was used with the neural network approximation. 

Method 3: This method replaced normal distribution with the Weibull function in method 2. The 
regression method was used for approximation. 

There was little difference between the mean value and deterministic value for stress by all three 
methods. The standard deviation was about 9 percent, which compared well with the 10-percent 
deviation in the primitive random variables. The design model had three random variables, consisting of 
one thickness along the web boundary, one in the plate central zone, and one in between (the 
intermediate region). To begin optimization calculations, the mean value of the variables were taken 
equal to the deterministic optimum solutions. The standard deviations were 10 percent of the mean value 
for all three of the design variables. CometBoards testbed was used for design optimization. Optimum 
weights for different probability of failure are given in Table 4. The rate ranged from a vulnerable 
design (with 999 failures in 1000 samples) to a reliable design with 1 failure in 1000 samples. 

Weight versus probability of success is plotted in Fig. 7 . There is some deviation in the weight 
calculated by the three different analysis methods as shown in the figure. For a 50-percent probability of 
success all three methods converged to 3.22 lbf for the optimum weight. For a reliable design with 1 
failure in 1000 samples the first method produced about a 10-percent higher weight than that obtained 
by the neural network technique. The difference reduced to 5 percent for the regression method. 
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Table 4. Optimum weight versus probability level for the web plate. 




Optimum weight of web plate 

(X) 

Probability 

Method 1 

Method 2 

Method 3 

(X failures in 

1000 

samples) 

level, p 

FPI and normal 
distribution 

Neural network 
and normal 
distribution 

Regression method 
and Weibull 
function 

(999) 

0.001 

1.16 

1.78 

1.75 

(975) 

0.025 

1.89 

2.26 

2.22 

(900) 

0.1 

2.33 

2.56 

2.54 

(700) 

0.3 

2.84 

2.94 

2.93 

(500) 

0.5 

3.22 

3.22 

3.22 

(300) 

0.7 

3.61 

3.52 

3.52 

(100) 

0.9 

4.21 

3.99 

3.98 

(25) 

0.975 

4.81 

4.47 

4.41 

(1) 

0.999 

5.96 

5.40 

5.16 
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Figure 7. Inverted-S graph for optimum weight versus probability level for web plate. 
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Figure 8. Basic reliability design concept; a qualitative illustration for induced stress, (a) Deterministic design, (b) 
Design for mean value, (c) Design for more than mean value, (d) Design for less than mean value. 
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8.2 THE INVERTED S-SHAPED GRAPH 


The weight versus probability of success for the web plate given in Table 4 is plotted in Fig. 7. 
Weight increased when reliability exceeded 50 percent. Weight decreased when the reliability was 
compromised. The weight versus reliability traced out the inverted-S-shaped graph. A design can be 
selected that depends on the level of risk acceptable to the situation. The illustration in Fig. 8 provides 
a simple explanation for the inverted-S shape of the graph. 

Consider a deterministic design calculated for an allowable stress limit of 25 ksi with a 100 lbf 
weight for the structure (see Fig. 8(a)). The structure is redesigned next for a 50-percent probability of 
success. Assuming a 1.5 safety factor in applied loads, the stress can be approximated to 17 ksi and the 
proportioned weight can be approximated as 67 lbf (see Fig. 8(b)). Let us consider a reliable design with 
1 failure in 100 000 samples. The stress value is likely to increase to 24 ksi because of an increase in the 
corresponding area of the distribution function (see Fig. 8(c)). The weight has to be increased to about 
80 lbf because of a high value for stress, as shown in Fig. 8(c). Consider next a failure-prone design with 
90 failures in 100 samples. The stress can be less, like 7 ksi, because of a reduced area under the 
distribution function (see Fig. 8(d)). The weight can be reduced to about 28 lbf because of the low stress 
level. In reliability-based design optimization, weight can become very heavy when reliability 
approaches unity; likewise, a lightweight design can be obtained when reliability is compromised. 





Figure 9. Model of raked wing tip of Boeing 767^100 extended range airliner. 
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8.3 RAKED WING TIP STRUCTURE 


Reliability-based optimization is illustrated considering the example of a composite airframe 
component of the Boeing 767-400 extended range airliner, the raked wing tip structure shown in Fig. 9. 
The structure is fabricated out of component parts made of metallic and composite materials. The 
available industrial design is referred to as the initial design. It is required to generate a deterministic 
optimum design as well as a reliability-based design for the component. The problem is treated in five 
separate sections: 

(1) Deterministic analysis 

(2) Deterministic optimization 

(3) and (4) Probabilistic analysis 

(5) Reliability-based design optimization 
8.3.1 DETERMINISTIC ANALYSIS 

The component has a wing-box type of construction with face sheets and web panels made of 
composites materials with an aluminum honeycomb core. Its components such as the root rib, tip rib, 
and leading edge parts are made of aluminum, and a set of titanium plates are also used. The 
MSC/Nastran code was used as the deterministic analysis tool. The finite element model had about 3000 
nodes and 17 000 degrees of freedom. The model was made of 3500 elements, consisting of CBAR, 
CBUSH, CELAS1, CHEXA, CONROD, CPENTA, CQUAD4, CROD, CSHEAR, and CTRIA3 
elements. For aluminum the Young’s modulus was in the range 10.2 to 10.9 million psi, the Poisson’s 
ratio was 0.33, and the density varied between 0.097 to 0.102 lbf/in 3 . Titanium had a Young’s modulus 
of 16 million psi, a Poisson’s ratio of 0.31, and a density of 0.16 lbf/in 3 . A typical Young’s modulus for 
composite laminate was 8 million psi with a 0.7 million psi for shear modulus and a density of about 
0.055 lbf/in 3 . A composite panel may contain between 6 and 40 plies, each of 0.009 in. thickness. The 
thickness for the aluminum honeycomb was in the range 0.5 to 0.8 in. with the shear modulus between 
2000 to 4000 psi and density of about 0.0017 lbf/in 3 . The structure was subjected to eight load cases. 


The CPU time for a single MSC/Nastran run was about 5 s for static analysis, and it required almost 1 
min to calculate a set of 20 frequencies on a Red Hat Enterprise Linux workstation, Release 4. Overall 
response of the initial model for the eight load cases is depicted in Table 5. 
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Table 5. Response of the composite raked wing tip. 


Load case 

Strain energy, 

u, 

in.-kip 

Maximum 

displacement, 

in. 

Maximum von 
Mises stress for 
metal, 

ksi 

Maximum principal 
microstrain in 
composite plies 

1 

9.3 

-5.63 

38.2 

2427 

2 

39.8 

11.56 

67.4 

5450 

3 

34.4 

10.59 

65.7 

4920 

4 

33.0 

10.36 

66.1 

4940 

5 

8.9 

-5.54 

35.5 

2600 

6 

54 

13.68 

87.7 

7470 

7 

44.3 

12.45 

84.8 

6980 

8 

28.2 

10.60 

61.6 

4180 

Frequency in hertz for first tl 

free modes: 12.69, 16.34 (cantilever mode], and 17.07 


The strain energy U in the first column in Table 5 was obtained as one-half of the product of load P 
and displacement X over all nodes using the following formula: 


u = 


J (strain x stress f/r = — 


Volume 


^ all nodes and directions' ^ 

I 

l M 


Z P J X J 


(18) 


The strain energies for the six load cases were in the range 9 to 54 in.-kip. Load case 6 contained the 
peak value for strain energy of 54 in.-kip. The maximum displacement for the load case was 13.68 in. 
Maximum strain in the composite components was 7470 microstrain. The maximum stress in metals was 
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87.7 ksi. The structure was designed for the critical load case 6 and checked against other load cases. 
The concept to design for a load case with maximum strain energy drastically reduced the number of 
calculations in design optimization, and it was subsequently proved correct for the problem. 

8.4 DETERMINISTIC OPTIMIZATION 

The objective of the optimization problem was to reduce the weight of the wingtip without changing its 
geometrical configuration. Only the thicknesses of the components were allowed to be adjusted within 
prescribed bounds. For design optimization the raked wingtip was separated into several substructures 
consisting of metallic and composite components. From a preliminary design optimization study it was 
concluded that the thickness of titanium components and aluminum parts like, root rib center section, 
leading and trailing edge webs, as well as leading edge and tip rib skins cannot be changed from a 
manufacturing consideration. Design variables associated with these parts were considered passive. The 
remaining four composite members were grouped to obtain a total of 13 active design variables. 
Grouping was on the basis of the number of plies. For example, the first design variable represented a 
proportional thickness parameter for the upper and lower skin panels. Initially the panel had 7 to 16 plies 
over honeycomb thicknesses ranging from 0.65 to 0.75 in. This variable contained a total of 211 
CQUAD4 and CTRIA3 elements. Likewise, the design variable 10 represented the area parameter for 
spars in the front web. The initial areas of the spars were in the range 0.18 to 0.393 in 2 and 92 CROD 
elements were used in the discretization. Other design variables were defined in a similar manner. 

For constraint formulation, the structure was separated into 203 groups of elements to obtain a total of 
203 strain constraints for the upper and lower skin panels as well as the spars. The composite rod 
elements were grouped to obtain 16 more strain constraints. Three translations and one rotation were 
constrained at a tip node of the structure for load case 6 as well as for load case 7 to obtain eight 
displacement constraints. The design model has a total of 227 behavior constraints. Constraint can be 
imposed on principal strain, for a failure theory, or on a strain component. The allowable strain was 
4000 microstrain. Displacement limitations along x-, y-, and z-directions were set at 0.5, 1.5, and 14 in., 
respectively. Maximum rotation was not to exceed 5°. Design optimization was performed using the 
testbed CometBoards. Deterministic optimum solution is given in Table 6. Only a normalized optimum 
solution can be given because of the proprietary nature of the data. 

Stress and strain obtained from the finite element model were presumed not to be accurate at some metal 
and composite interface locations. Such localized regions were avoided in design calculations by 
excluding the associated strain constraints. In Table 6 column 2, model 1, refers to a configuration that 
was obtained by excluding the strain constraints for a set of eight interface elements. Likewise model 2 
excluded strain constraints associated with a total of 24 elements. Generation of an optimum solution for 
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model 1, for example required about 39 CPU minutes. Optimum weight for model 1 was 16 percent 
lighter than that of the original design. The weight saving was 20 percent for model 2. Design changed 
throughout the structure except for variables 12 and 13, which represented the minimum thicknesses. 
Maximum reduction was observed for design variable 8 with a 44 percent change. The thickness 
reduced to 0.44 in. if the original value was 1 in. In other words, in the original structure the entire 
region that was associated with design variable 8 represented an overdesigned condition. The thickness 
variable 8 for model 1 increased to 1.24 in. from its assumed value of 1 in., or the region originally was 
underdesigned. The design process redistributed the strain field with many active constraints. There was 
little change in the displacement state or frequency value for model 1. For model 2, maximum 
displacement reduced by 2.5 percent, but its frequency increased by 25 percent because of a 20 percent 
reduction in its weight. 
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8.5 PROBABILISTIC ANALYSIS 

Probabilistic analysis requires a geometrical model, a load model and a material model: 


Geometrical model: All 13 design variables of the deterministic design optimization were considered to 
be random variables. Their mean values were equal to the deterministic optimum design solutions. Their 
standard deviations were set at 7.5 percent of the mean values. The thicknesses of the honeycomb were 
considered as random variables. Their mean values were equal to that of the initial design with a 
standard deviation of ten percent of the mean values. Cross-sectional areas of the bars were considered 
to be deterministic parameters and were equal to the optimum solutions. 


Table 6. Normalized deterministic optimum solution 


Design variables 

Weight 

Original design 

Design model 1 

Design model 2 

100 units* 

84.0 units 

80.0 units 

Normalized 
to unity 

Change in percent 
(100 represents no change) 

1 

1.0 

70.1 

70.1 

2 

1.0 

104.5 

99.83 

3 

1.0 

110.3 

111.83 

4 

1.0 

75.64 

76.97 

5 

1.0 

44.85 

44.45 

6 

1.0 

123.35 

80.96 

7 

1.0 

84.17 

94.45 

8 

1.0 

44.08 

44.88 

9 

1.0 

86.12 

83.59 

10 

1.0 

88.09 

88.09 

11 

1.0 

73.83 

73.83 

12 

1.0 

100.00 

100.00 

13 

1.0 

100.00 

100.00 

Active strain 
constraints 

— 

6 with 4000 
micro strain; 9 
exceeded 3600 
micro strain 

6 with 4000 

microstrain; 9 exceeded 
3600 micro strain 

Displacement in 
z-direction, in. 
(limit was 14 in.) 

13.68 

13.14 

13.35 

Rotation, deg(limit 
was 5°) 

4.24 

4.48 

4.72 

Frequency, Hz 

16.36 

16.45 

20.33 


^Normalized to 100 units. 
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Load model: Each load component is considered to be a random variable with a mean value equal to 
66.67 percent of the deterministic value, with a 10-percent standard deviation. For example, if a 
deterministic load component is 100 lbf, then its random counter part would have 66.67 lbf for its mean 
value with 6.67 lbf for the standard deviation. Each load component was changed in a similar manner. 
This refers to load model A. A second load model was recommended and added. Thus, there were two 
load models, load model A and load model B. Load model B was generated following the procedure for 
load model A but with a reduced mean value of 50 percent of the deterministic load value. A standard 
deviation was set at 10 percent of the mean load value. 

Material model: Modulus of elasticity as well as Poisson’s ratio were considered to be random 
variables. The standard deviations for all the material parameters were set at 7.5 percent of their mean 
values. The mean values were set to 105 percent of their deterministic values for the Young’s modulus and 
shear modulus of fabric as well as the shear modulus for the honeycomb. The mean value for Poisson’s 
ratio was set to 100 percent of its deterministic value. 

8 . 5.1 Probabilistic Solution for Strain and Displacement 

The probabilistic solution was obtained using the MSC/Nastran code and the FPI module of the 
NESSUS software. Four different types of distributions were considered: normal distribution, Weibull 
function, lognormal, and Gumbel 22 type 1 distributions. The probabilistic analysis produced the mean 
value of strain (in the element number 11531) to be 5050 microstrain by all four types of distribution 
functions (see row 1 in Table 7). Likewise, the standard deviation remained the same at 16 percent of the 
mean value by all four distribution functions. The mean value of displacement (at node 1 1243 in the z- 
direction) was 9.6 inch with a standard deviation of 10 percent by all four distributions. The observation 
that the mean values and standard deviations did not change for different distribution types was along 
the expected line. 

Probabilistic solution for the strain and the displacement for four different types of distributions for 
different failure rates is given in Table 7. For a 50-percent probability of failure, or one failure in two 
samples, the mean value and the calculated strain are equal at 5050 microstrain for normal distribution. 
The same is not true for Weibull, lognormal, or Gumbel type 1 distributions because of a lack of 
symmetry in such functions. Consider next, one failure in one million samples. The Weibull prediction 
was conservative. It was about 85 percent of the normal distribution for strain as well as displacement. 
Estimates by Gumbel type 1 distribution was on the higher side. For one failure in a million samples, it 
predicted about 40 to 45 percent higher strain and displacement than that for the normal distribution 
function. 


NASA/CR— 201 3-2 1 7748 


198 



Table 7. Probabilistic response for different probability levels.* 


N samples 

Microstrain in element 11531 

Displacement at node 1 1243 
along z-direction 


Normal 

Weibull 

Lognormal 

Gumbel 
type 1 

Normal 

Weibull 

Lognormal 

Gumbel 
type 1 

Mean 

value 

5050 

5050 

5050 

5050 

9.7 

9.7 

9.7 

9.7 

Standard 

deviation 

810 

810 

810 

810 

0.98 

0.98 

0.98 

0.98 

2 

5050 

5123 

4986 

4917 

9.7 

9.7 

9.6 

9.4 

1 000 

7555 

6967 

8162 

9050 

12.6 

11.7 

13.1 

14.4 

50 000 

8379 

7398 

9599 

11523 

13.6 

12.2 

14.5 

17.4 

100 000 

8507 

7459 

9843 

11961 

13.7 

12.2 

14.7 

17.9 

500 000 

8903 

7591 

10402 

12977 

14.1 

12.3 

15.2 

19.1 

1 million 

9264 

7644 

10641 

13415 

14.2 

12.4 

15.5 

19.7 

10 

million 

9264 

7804 

11424 

14870 

14.7 

12.6 

16.2 

21.4 


The first and second rows provide the mean values and standard deviations. Remaining 
data in Table 7 represent mean values. 
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8.6 PROBABILISTIC DESIGN OPTIMIZA TION 

Probabilistic optimization is performed for the two deterministic design models, model 1 and model 2 
as discussed earlier (see Table 6). The probabilistic design calculation used the information given for 
stochastic analysis along with additional data required to formulate failure modes or the design 
constraints. It was assumed that mean value of allowable strain was 25percent higher than its 
deterministic limit, while the standard deviation was set at 7.5 percent of the mean value. For example, 
the limits for strain were set with a mean value of 5000 microstrain and a standard deviation of 375. 
From deterministic optimization calculations, it was observed that only the displacement limitation 
along the z-direction and rotation had some influence in the design while other stiffness constraints were 
quite passive. Thus, for probabilistic design only two stiffness constraints were retained. At node 11243, 
along the z-direction, the mean value and standard deviation for the displacement limit were 17. 1 and 
1.71 in., respectively. The mean value and standard deviation for rotation limit were set at 6.25° and 
0.625°, respectively. For the design calculation normal distributions were assumed for all random 
parameters. 

The optimization calculation required continuous running of the code for more than 5 days, but the 
execution was smooth and eventless. The stochastic optimum solution was quite similar to that obtained 
for deterministic optimization depicted in Table 6. There were numerous active strain, displacement, and 
rotation constraints. Normalized optimum weight obtained for different failure rates is given in Table 8. 
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Figure 10. Half of inverted-S graph for design model 1 and load model A. 
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Asymptotic to infinity 


Table 8. Probability of failure versus weight 


N 

samples 


Normalized optimum weight 


Design model 1 and 
load model A 

(strength + 
stiffness) 

Design model 1 and 
load model A 

(strength only) 

Design model 1 and 
load model B 

[strength or 

(strength + 
stiffness)] 

Design model 2 
and load model A 

(strength + 
stiffness) 

2 

64.68 

64.68 

62.24 

63.20 

10 

67.43 

67.43 

63.57 

64.36 

100 

70.47 

70.47 

65.51 

66.97 

1000 

73.75 

73.75 

66.97 

69.57 

10 000 

76.64 

76.64 

68.23 

72.48 

100 000 

79.28 

79.03 

69.70 

76.25 

1000 000 

92.93 

82.27 

71.57 

91.21 

1250 000 

94.88 

82.55 

71.76 

93.41 

2 000 000 

100.00 

83.15 

72.16 

98.60 
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The normalized optimum weight was set to 100 for strength and stiffness constraints for load case A 
in design model 1 for 1 failure in 2 million samples. This design exhibited nine active constraints 
consisting of eight strain and one displacement limitation. For an active stochastic constraint the first 
factor in Eq.(7) (based on mean values) was equal to the second factor (which was a function of standard 
deviation and the phi critical parameter) in magnitude but with opposite sign. The frequency was 15.94 
Hz for the optimum solution. One stochastic design optimization run with 61 p levels required 128 h, or 
5.33 days. 


Both stress and 



1 10 100 1000 10 000 10 00001 000 000 10 000 000 (log N) 

N samples 


Figure 11. Inverted-S graph in log scale for design model 1 and load model A. 


8.7 THE INVERTED S-SHAPED GRAPH 

The weight versus probability of success given in Table 8 is plotted in Figs. 10 and 11. In Fig. 10, the 
x-axis represents N (as in one failure in N samples) and it begins at N = 2, or 50 percent probability of 
success. This graph represents one-half of the inverted-S graph because probability less than 50 percent 
is not included. This graph is for load case A, design model 1 ; both strength and stiffness constraints are 
considered, and column 1 in Table 8 contains the weight information. 
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The same information is replotted in Fig. 11, but a logarithmic scale, or log(W), is used along the x- 
axis. This figure can be approximated into two linear segments. At the intersection point (SD) both 
strength and stiffness constraints are active. From the origin to the point SD, only strength constraints 
were active. From point SD onwards, both strain and stiffness constraints are active. 

The variation of weight shown in Table 8 was as expected. The weight is least for a 50-percent rate of 
success, or for N= 2. The weight increased when failure rate was reduced. The weight in the first two 
colu mn s in the Table 8 coincided up to 1 failure in 10 000 samples because only strength constraints 
were active. The weight increased when both strength and stiffness constraints became active (see 
column 2, in Table 8). Weight shown in column 4 was lighter than that in colu mn 2 because load 
model B was less severe than load model A. Weight shown in column 5 was lighter than that in column 
2 because design model 2 was obtained by removing few severe strain constraints from model 1. 



Deterministic design Stochastic design 

Industrial design: optimization: optimization 

Weight = 100 units Weight = 84.5 for one failure in 

2 million samples: 

Weight = 83 

Figure 12. Optimization process more evenly distributed the strain field in wing tip. 

The distribution of the strain state in the wing tip is illustrated in Fig. 12 for three cases. The first case 
is the initial design that was obtained by the industry. The second one is the deterministic solution while 
the third one represents the stochastic design. The strain distribution is uneven for the initial design, and 
strain exceeds 4000 microstrain at a few localized locations. The distribution of the strain state is more 
even for the deterministic as well as the stochastic design solutions. 
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Legend 

Case 1: Original design 

Case 2: Deterministic optimum design for design model 1 and load model A 
Case 3: Deterministic optimum design for design model 2 and load model B 
Case 4 and 5: Stochastic design for design model 1 and load model A 
Case 6: Stochastic design for design model 2 and load model B 


84.5 


— Standard deviation 




1.26 


94 



91.7 


Strain and Strain and (1 in 1 000 000) (1 in 1 000 000) (1 in 1 000 000) 
displacement displacement Strain only Strain and Strain and 

displacement displacement 

Design case 1 2 3 4 5 6 


Figure 13. Optimum weight for six design cases. 


Optimum weight for six design cases are depicted in Fig. 13 in a bar chart. The first case with a 
weight assumed at 100 units represents the original design. All five calculated design cases have lower 
weight than the original design. The deterministic optimum solution (for strength and stiffness 
constraints with design model 1 and load case A) generated a 15.5-percent lighter design than the 
original design. The weight decreased further to 80.9 units when the stiffness constraint was relaxed, see 
design case 3. Stochastic optimization produced a mean value for the weightof 83 units with a standard 
deviation of 1.3 units for strength constraints for design model 1 and load case A for 1 failure in 1 
million samples, (see design case 4). The mean value for the weight and a standard deviation increased 
to 94 and 1.26 units when stiffness constraint was added to strength limitation, (see design case 5). The 
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mean value of weight and standard deviation reduced to 91.7 and 1.23 units, respectively, for strength 
and stiffness constraints for design model 2 and load case B (see design case 6 in Fig. 13). Over all, 
weight can be reduced up to 20 percent depending on the choice for model for design, load, and rate of 
failure. The standard deviation was small and remained at about 1 percent for the three design cases 4 to 
6. 


Sensitivity analysis was performed for the principal strain (in element 11531) for deterministic as 
well as stochastic design optimization methods. The element strain was most sensitive to load and elastic 
modulus E for fabric as well as for design variable 6 (see Table 6) that contained element 11531. The 
strain was not sensitive to the other 19 random variables. Both deterministic as well as stochastic 
methods identified the same set of variables, namely DV6, load, and Young’s modulus. 


Table 9. CPU time for design and analysis of wing tip. 


Activity 

CPU time 

One static analysis cycle in MSC/Nastran 

5 seconds 

Dynamic analysis to calculate 20 frequencies in 
MSC/Nastran 

51 

seconds 

Deterministic optimization for design model 1 

39 

minutes 

Stochastic analysis 

47 

minutes 

Stochastic optimization: Design model 1 and load 
model A 

128 

hours 

Stochastic optimization: Design model 2 and load 
model A 

126 

hours 
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The CPU time to solution is given in Table 9. The calculation used CometBoards, which is NASA in 
house software, along with MSC/Nastran version 2005.5.0 (2005R3) and the FPI module of NESSUS 
level 6.2 code (1995). Calculations used a Red Hat Linux 2.6.9-67.ELsmp O/S, with x86_64 
architecture, 2600 MHz, 4 cpu, 8 GB of memory, and 32-bit numeric format. 

One static analysis cycle required about 5 CPU seconds. The run time increased to 51 s for dynamic 
analysis. Stochastic analysis required about 47 min. Deterministic optimization required about 39 min. 
The CPU time for stochastic optimization was enormous at 126 to 128 h of continuous calculations. 

9.0 Concluding Remarks 

For finite element structural analysis in the probabilistic domain, analytic expressions have been 
derived for the mean value and covariance matrix of forces and displacements to quadratic order. 
Expressions have been derived for both the Integrated Force Method (or flexibility method) and the 
Dual Integrated Force Method (or stiffness method) using a perturbation technique. The derivation 
included mechanical and initial deformation loads, variation in material properties as well as design 
parameters. The mean of the force and displacement contain zeroth order and quadratic terms. The 
zeroth, linear and quadratic terms of expansions contributed to the covariance for both variables. The 
expressions have been programmed in closed form in Macsyma for simple structures. For the structures 
subjected to thermo-mechanical loads, the stochastic mean changed little, but the covariance itself was 
significant. The force response was most sensitive to the mechanical load and least sensitive to Young’s 
modulus. The displacement response is most sensitive to the Young’s modulus and a mechanical load 
component. 

The test bed CometBoards with MSC/Nastran and a fast probability integrator successfully generated 
reliability-based design optimization for an industrial strength raked wing tip structure of a Boeing 767- 
400 extended range airliner made of metallic and composite materials. The optimization run that 
required 128 h or 5.25 days of continuous execution was eventless. The optimization methodology 
requires probabilistic models for load, material properties, failure theory, and design parameters. 
Accuracy of the design solution depends on the models. Stochastic optimum weight versus reliability 
traced out an inverted-S-shaped graph. Weight increased when risk was reduced and vice versa. The 
inverted-S graph degenerated into linear segments when a logarithmic scale was used for the x-axis. The 
optimization process redistributed the strain field in the structure and achieved up to a 20-percent 
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reduction in weight over traditional design. Optimum weight was comparable between a deterministic 
solution and a highly reliable stochastic design. Inclusion of standard deviation as an additional 
objective function should be examined further because it offers very little variation and has insignificant 
impact. It may be a manufacturing issue with little relevance to analytical design calculations. Design of 
an aircraft structural component has to be obtained for multiple load conditions. The maximum strain 
energy criteria identified a few critical load cases from the many load conditions. The design generated 
for the critical load was satisfactory. This approach reduced calculations to a very great extent. The 
design sensitivity can identify critical zones for redesign consideration. Both deterministic and 
stochastic concepts identified identical zones. There was no preference to either concept. 
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